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We review the state of the art of the theory of Euclidean random matrices, focusing on the density 
of their eigenvalues. Both Hermitian and non-Hermitian matrices are considered and links with 
simpler, standard random matrix ensembles are established. We discuss applications of Euclidean 
random matrices to contemporary problems in condensed matter physics, optics, and quantum 
chaos. 
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I. INTRODUCTION 

Random matrix theory (RMT) is a powerful tool of 
modern theoretical physics ( Mehtal 2004 1. Its main goal 



is to calculate the statistical properties of eigenvalues 
and eigenvectors of large matrices with random elements. 
First introduced by Wishart (1928) and then used by 



Wigner ( 1955 ) to describe the statistics of energy levels 



in complex nuclei, random matrices are nowadays om- 



nipresent in physics ( 


Beenakker 


1997 


Brody et al. 


1981 


Guhr et al. 


1998 


Tu 


ino and Verdii 


2004) 


. The majority 


of works — including the seminal papers by 


Wigner 


1955) 



-deal with Hermitian matrices. 
Hermitian matrices are of special importance in physics 



2 



because of the Hermiticity of operators associated with 
observables in quantum mechanics. On the other hand, 
non-Hermitian random matrices have also attracted con- 



siderable attention (Feinberg 20061 iFeinberg and Zee 



1997a b| |1999a|b||Hatano an d Nelson||1996||Janik et al. 



1997a b; Jarosz and Nowak 2006). They can be used to 



model such physical phenomena as scattering in dissipa- 



tive or open systems ( 


Fyodorov et al.\ 


2005 


Fyodorov 


and Sommers 


1997 12003 


Haake et al.\ 


1992), dynamics 


of neural networks (|Rajan and Abbott 


2006 


Sommers 



et al. 1988), diffusion in random velocity fields (Chalker 



and Wang 1997), or chiral symmetry breaking of the 



QCD Dirac op erator ( Stephanov 1996 Verbaarschot and 
Wettig] 12000). 



A special class of random matrices are the so-called Eu- 



clidean random matrices (ERMs) (Mezard et al. 1999). 
The elements A^ of a N x N Euclidean random ma- 
trix A are given by a deterministic function / of posi- 
tions of pairs of points that are randomly distributed in 



a finite region V of Euclidean space: A 



*7 



i = 1,...,N. Hermitian ERM models play an impor- 
tant role in the theoretical description of vibrations in 



topologically disordered systems (Ganter and Schirma- 



cherH2011; Grigera et a/.||2011||2001a||2002H2003||2001b 


Mezard et al. 


1999 1 and relaxation in glasses ( Amir et al. 


2008 


2010 


2011 


). They have been used to study Ander- 



son localization (Bogomolny et al. 2003 Ciliberti et al. 



2005 Huang and Wu, 2009). A number of analytic ap- 



proaches were developed to deal with Hermitian ERMs 



( Amir et al. 


2010| Bogomolny et al. 2003 Chamon and 


Mudry 2001 


Ciliberti et al. 2005| Ganter and Schirma- 


che4|2011; Grigera et aZ.||2011| 2001a, 2002, 2003, 2001b| 


Mezard et al 


[ |1999 Skipetrov and Goetschy[|2011). In 



contrast, an analytic approach to non-Hermitian ERMs 



was pro posed only recently (Goetschy and Skipetrov 
2011a|b |. A particular non-Hermitian ERM — the so- 
called Green's matrix — was a subject of extensive nu- 



merical studies (Antezza et al. 2010| Pinheiro et al. 



2004||Rusek et aL||2000a[|Skipetrov and Goetschy[|201ip. 



The principal difficulties that one encounters when trying 
to develop a theory for ERMs stem from the nontrivial 
statistics of their elements and the correlations between 
them. Both are not known analytically and are often dif- 
ficult to calculate. This is in sharp contrast with standard 



2004 ) in which the joint probability distribution of the el- 
ements of the random matrix under study is the starting 
point of analysis. 

The purpose of this paper is to review the existing an- 
alytical results for the eigenvalue density of large ERMs 
(N 3> 1) and to discuss their applications to modern 
problems in condensed matter physics, optics, and quan- 
tum chaos. Analytical results will be compared with 
numerical diagonalization and special attention will be 
given to the range of validity of the former. To derive 
the analytical results, we will use both the diagrammatic 
techniques which are more familiar to a theoretical physi- 
cist, and the recently introduced free probability theory 



(Voiculescu et al. 1992[ ) which is a mathematically in- 
volved but powerful technique. Arriving at the same re- 
sult in two different ways will help to understand which 
approximations were essential and which were not and 
will, we hope, promote the use of the new mathemati- 
cal tool — free probability theory — by the community of 
physicists. 

Because ERMs are a relatively new class of random 
matrices, we believe that it is essential to make links 
with more 'standard' random matrix ensembles, such as 
the matrices with independent, identically distributed 

We will 



2004) 



elements or Wishart matrices (Mehta 
show that for some of the ERMs under study and in 
a restricted range of parameters, the eigenvalue density 
is very close to that of the known, simpler matrix en- 



sembles: Wigner semi-circle (Mehta 2004), Marchenko- 
Pastur law (Marchenko and Pastur 1967), Ginibre's disk 



(Ginibre 1965 Girko 1985). These well-known distribu- 
tions will be introduced and discussed as well, to ensure 
a self-contained presentation. 

Finally, we will discuss what is known about the prop- 
erties of eigenvectors of Euclidean random matrices. The 
study of eigenvectors is a much more difficult problem, 
which explains our very poor understanding of their sta- 
tistical properties. We will focus on the most interesting 
question of spatial (Anderson) localization of eigenvec- 
tors. 

The paper contains three main sections: Sec. [il] is de- 
voted to Hermitian matrices, Sec. |III| — to non-Hermitian 
ERMs, and, finally, Sec. |IV| discusses applications of 
ERMs in physics. Sections |II| and |TTT] are self-contained, 
even though the latter makes frequent references to the 
former to help the reader understand the more compli- 
cated case of non-Hermitian matrices by analogy with 
the simpler Hermitian one. This review is supposed to be 
accessible to everyone with a general theoretical physics 
background and does not require knowledge of RMT to 
be understood. It can serve as an introduction to RMT, 
even though the view it provides is biased by its main fo- 
cus on ERM ensembles which are only a small fraction of 
many important ensembles of RMT. For a more general 
and unbiased view of the modern RMT, we invite the 
interested reader to consult the book of iMehtal (12004) or 



approaches ( 


Beenakker 


1997 


Guhr et al. 


1998 


Mehta 


et al. ( 


2001) 



RMT and of their applications in physics can be found in 



review papers by Beenakker (1997); Brody et al. (1981); 
Mitchell et al.\ 120101). 



II. HERMITIAN EUCLIDEAN RANDOM MATRIX 
THEORY 

In this section, we start by introducing a number 
of well-known ensembles of Hermitian random matrices 



(Sec. II. A), including the Euclidean random matrix en- 
semble that will be of interest for us (Sec. |II.A.3 ). We 
then review the main approaches to study the statistical 



3 



properties of eigenvalues of Hermitian random matrices 
in general and of Euclidean matrices in particular (Sees. 
II.BHII.HI) and apply them to analyze the eigenvalue dis- 



tributions of several specific Euclidean random matrices 
(Sec.|nj|. 



A. Some standard random matrix ensembles 

1. Gaussian matrices 

The best known random matrix ensembles are proba- 
bly the Gaussian ensembles. They are ensembles of iVx N 
Hermitian matrices A = , that have independent and 
identically distributed (i.i.d.) zero-mean Gaussian en- 
tries. The probability distribution of A is 



P(A) = C N e~- 



■TtA 



(1) 



where Cn is a normalization constant, and f3 is the sym- 
metry index, that counts the number of degrees of free- 
dom in the matrix elements. 

For our purpose, it is sufficient to consider matrices A 
with entries being either real or complex numbers (0 = 1 
or 2). Let us first analyze the case of complex elements, 
for which (3 = 2. Since the transformation A — > UAU -1 , 
with U unitary, leaves P(A) invariant, the ensemble is 
called 'Gaussian unitary ensemble' (GUE). From Eq. (JlJ, 
we easily verify that the second moments of Ay are 



(AijAki) 



GUE 09 = 2). 



(2) 



On the other hand, if elements of A are real num- 
bers, (3 = 1, and the transformation A —> UAU^ 1 leaves 
P(A) invariant for U orthogonal. The ensemble is called 
'Gaussian orthogonal ensemble' (GOE), and the second 
moments are given by 

(AijAu) = ^(daSjk + 6 ik 5ji), GOE (ft = 1). (3) 

As we shall see later, the density of eigenvalues of a 
Gaussian matrix A converges, in the limit N — > oo, to 
the so-called 'semicircle' law, first discovered by |Wigner| 
p55l). 



2. Wishart matrices 



Another ensemble of particular interest for us is the 



Wishart ensemble, that is as old as RMT itself (Wishart 



1928). It is useful in many contexts, such as neural 



networks, image processing, or wireless communications, 
where Wishart matrices naturally arise to characterize 
the singular values of 'channel matrices' rtTulino and 



Verdu 



2004). A N x N Wishart matrix A is of the form 



A 



(4) 



where H is a rectangular N x M matrix, with columns 
that are zero-mean independent real or complex Gaus- 



and Verdu 2004). In this section, we will work with H 
complex and X proportional to the identity matrix In- 
In this case, the entries of H are zero-mean i.i.d. complex 
Gaussian random numbers. The probability distribution 
of the non-Hermitian matrix H is 



P(H) = C N , M e 



-NTiHH f 



so that the second-moments of elements of H obey 



(H ia H 



Pi' 



N 



ijo a p - (H ai Hjp) 



(5) 



(6) 



For c = N/M < 1, |Wishart| ( |1928[ ) showed that the 
probability distribution of Q is given by 



P(A) = C' NtM detA M - N »- NTrA 



(7) 



see also ( |Tulino and Verdu" 2004). Quite surprisingly, no 
such explicit formula was known for c > 1 ('anti- Wishart 
case') until recently ( Janik and Nowak |2003 ). However, 
as far as the eigenvalue distribution of A is concerned, it 
is straightforward to obtain the result for c > 1 from the 
one for c < 1 (see Sec. 



In Sec. |II.C.4[ we wi 



II.C.4) 



1 see that the eigenvalue distribu- 
tion of A = HH^ converges, in the limit N, M — > oo with 
c = N/M fixed, to the so-called Marchenko-Pastur law. 
It was first established by Marchenko and Pastur ( 19671, 



and then rediscovered several times (Tulino and Verdu 



2004). 



3. Euclidean random matrices 

As explained in the introduction, ERMs are matrices 
with elements Ay defined with the help of some deter- 
ministic function / of positions of pairs of points: 



Ai 



(8) 



where we introduced an operator A associated with the 
matrix A. We assume that the N points rj are randomly 
distributed inside some region V of d-dimcnsional space 
with a uniform density p = N/V . In this review, we 
will mainly focus on d = 3. Most of the analysis that 
will be presented below can be extended to the case of 
correlated positions of points by replacing /(r^iv,-) with 
f(rj, rjWr|— rA where ff(r) is the pair correlation func- 



tion ( Martin-Mayor et all 



2001) 



Contrary to Gaussian or Wishart matrices, the proba- 
bility distribution P(A) is not known analytically. When 
computing the statistical properties of A, averaging is 
performed not with respect to P(A) but with respect to 
the probability density of points {r^} in V . 

Depending on a particular physical problem under 
study, the matrix A defined by Eq. ^ may be subject 
to additional constraints. For ERMs that arise in the 
study of vibrational modes of an amorphous solid, in- 
stantaneous normal modes of a liquid, or random master 



sian random vectors with covariance matrix X (Tulino equations (Mezard et al. 1999), a condition Ay = 



4 



expresses the global translational invariance (conserva- 
tion of momentum in the case of propagating excita- 
tions). Typically, it is obeyed by using Eq. Q only for 
off-diagonal elements of A and adopting a different defini- 
tion for diagonal elements, An = — Ylj^i Aij- Both cases 



can be combined in a single definition (Mezard et al. 



1999) 



Ai 



f{ r i,*j) 



N 

u5 i3 ^2f{Y u Y k ) 

k=l 



(9) 



where u — or 1. 

A useful trick to study statistical properties of matrices 
defined by Eq. ^ is to chang e the basis from {r^} to 
hj)a}, which is orthogonal in V ( Skipetrov and Goetschy| 
|201l|). Inserting the closure relation 1 = ^ Q \')pa}('4'a\ 
into Eq. d8|), we obtain for arbitrary V: 



where 



A = HTH\ 



Hia = —pz{ri\^a), 

T a /3 = p(ip a \A\ip0}- 



(10) 

(11) 
(12) 



In Eq. (11) and (U2|, the prefactor p is introduced for 



later convenience. In a rectangular box, for example, 
\ipa) = |k Q ) with (r|k Q ) = exp(ik a r) / \/V , so that T a p 
are simply the Fourier coefficients of /(rj, rj): 



'-a/3 



N 



d d r t d d r, 



/( r i,ij)e 



-i(k a 



- k /3Tj) 



(13) 



The advantage of the representation ( 10 ) lies in the sepa- 



ration of two different sources of complexity: the matrix 
H is random but independent of the function /, whereas 
the matrix T depends on / but is not random. 
Furthermore, we assume that 



/ d d r </> a (r) = 0, 
Jv 



(14) 



which in a box is obeyed for all a except when k Q = 
0. We readily find that Hi a are identically distributed 
random variables with zero mean and variance equal to 
1/N: 



(Hi a ) 

(H ia H* jj3 ) 



1pa(*i) = 0, 



(H ia H 



1 f d a n 

VpJv v 

1 ff d d r, d rf r, 
pJJv V V 

(H ia ){H*p) = (i^j), 

1 r d d r, I 



i/3/ 



pJv V 



1p a (Ti)lp}(Ti) = —S af s 



(15) 

(16) 
(17) 



Equations (16) and (17) show that H satisfies Eq. ([6]), 
i.e. that the covariance matrix of the columns of H is 



S = Iff/N. If Hi a were Gaussian random variables, 
then this would be sufficient to conclude that Hi a are 
independent. However, they are not Gaussian and hence 
not necessarily independent. For example, the cumu- 
lant (AijAjiAij) c is not zero. It turns out that neglect- 
ing these complications and assuming Hi a Gaussian i.i.d. 
amounts to disregarding a class of 'dependent scatter- 
ing' events corresponding to the formation of 'cavities' 
by clusters of points rj. 

In Sec. |II.F| we will explicitly assume that Hi a are in- 
dependent Gaussian random variables. This assumption 
largely simplifies calculations but may limit the applica- 
bility of results, at least for certain types of Euclidean 
matrices, as we will see later. Within this assumption, 
the only but crucial difference between an ERM ( 10 1 and 
a Wishart matrix Q is the matrix T that contains all 
information about the function / defining the ERM. It 
can modify the eigenvalue distribution in a non-trivial 
way and lead to transitions between topologically differ- 
ent supports V of the eigenvalue density. Illustrations of 
such transitions are given by the examples considered in 
Sees. HTL3] and |HZ1 



B. Resolvent, Blue function, and 7£-transform 

Eigenvalues A„ of a N x N Hermitian matrix A are 
real. Their density, 



1 / N 

KA) = ^(£<Ka-a„) 



(18) 



\n=l 

can be obtained from the (one-point) resolvent 

-It— ' 



g(z) = i ( Tr^-p 
aw N \ z- A 



N \^z-A r , 

\n— 1 



(19) 



Using the standard relation lim e _ > . + l/(A + ie) = Pl/A — 
iwS( A) ( P denotes the Cauchy principal value), we rewrite 
Eq. flM as 



g(A + ie) = P J 



dA'l^-iTr^A), (20) 



so that p(A) may be reconstructed from cither the imag- 
inary or the real part of <?(A + ie): 1 



P (A) 



— lim InWA + ie), 

7T e-j-0+ 



I' /°°dA'|^=Re 5 (A + ie). 



(21) 
(22) 



1 Physically, g(A + ie) is the Fourier transform of the causal prop- 
agator e~ lAt S(t) [with 0(i) the Heaviside step function], and 
therefore its real and imaginary parts obey Kramers-Kronig re- 
lations. 



5 



Useless to say, Eq. (21 ) provides a much more direct way Both of them are fundamental objects of the free random 



to compute p(A) than the Fredholm integral equation of 
the first kind (22 1. However, the inversion of the lat- 



ter may sometimes yield the solution in a very efficient 
manner. Indeed, if p(A) has a finite support [a, b], the so- 



lution of Eq. (22 1 is given by Tricomi's theorem ( Tricomi 
1957}: 



p(A) 



(23) 



7rV(A-a)(6-A) 

" v /(A^- a )(6-AQ p 
,t P / dA — _ Reg (A + le) 



Such an expression for p(A) turns out to be particularly 
useful within the framework of the Dyson gas model (see 
SecllLCb. 



In order to compute g(z), we can rewrite Eq. (19) in 



various forms. Each of them is the starting point of a spe- 
cific analysis developed in the following sections. First, 
we note that 



N 

E 



i 



A. 



= d, In 



N 



n^- A «) 



and express the resolvent as 



1 



j(z) = -d, (In det(z - A)) 



(24) 



(25) 



This expression will be used in the field-theoretical ap- 
proach presented in Sec. |II.E) Another interesting expres- 
sion for g(z) is a decomposition in terms of the moments 
ofp(A), 



(A") = J°° dAp(A)A" = i<TWn 

J — OO ^ * 



(26) 



For Hermitian matrices, g{z) is a holomorphic function 
of z G C except for some cuts along the real axis where 
eigenvalues of A are concentrated. Therefore, we can 
reconstruct g(z) for all z by analytic continuation of its 
series expansion 



^ (A") 

#( z ) = 



(27) 



n=0 



which, in general, converges only in the vicinity of 



variable theory, discussed in Sec. II. G In particular, 72(z) 



is the generating function of the 'free cumulants' (see Sec. 
II. G for more details). According to Eq. (28), B{z) and 
TZ(z) are related to the self-energy a(z) by 



a(z)=n[g(z)}, 
B{z) = - + a[B(z)] 



(31) 
(32) 



Let us now mention a couple of properties useful for 
further analysis. The functions g(z), B(z), and 72(z) obey 
the following scaling relations: 



g a A{z) 

B a A{z) 
TZ aA (z) 



1 



9A(z/a), 



oSa{o-z\ 
aR-A{az), 



(33) 



where a e C*. Besides, the moments (A") can be ob- 
tained from g(z), B(z), and TZ(z). Using Eqs. ([27}, ([29}, 
and ( 30 ) , we readily show that 



(A") 
(A") 



d n+1 g(z) 



(n + 1)! d(l/z)™+ 1 
B 2 (z) d 



1 



(n + 1)! 
In particular, 



B'(z) dz 



B 2 (z) 
B'{z) 



(A)=^(0), 
varA=((A-(A)) 2 )=ft'(z)Uo, 



(34) 
(35) 



(36) 
(37) 



where B'{z) = dB(z)/dz and Tl'(z) = dTZ(z)/dz. Fi- 
nally, we note that the boundaries A* of the domain of 
existence of eigenvalues, p(A») = , are given by the fol- 
lowing simple relations (Zee 1996): 



<?'(A*) - 

S'(A*) : 

C. Mapping to the Dyson gas 



> oo, 
0. 



(38) 
(39) 



oo. We will work with the representation (27) in Sec. II. F 1. Dyson gas picture 



to perform a diagrammatic computation of g(z). In this 
perspective, it is also convenient to define the self-energy 
rj(z), that contains all irreducible diagrams in Eq. p7|: 



9(z) 



1 



a(zY 



(28) 



Observing that the electric field created by a point 
charge in two dimensions is inversely proportional to the 
distance from the charge, we can interpret the resolvent 



Other important objects for us are the functional in- 
verse of g{z), also called the Blue function, and the 72.- 
transform: 



B(z) 



72(z) = B{z) 



(29) 
(30) 



( 19 1 as the electric field created at a point z in the com- 
plex plane by charges q = 1 situated at positions A„ on 
the real axis. This suggests an analogy between the sta- 
tistical properties of random matrices and those of a gas 
of charged particles rest ricted to move in one dimension , 
the so-called Dyson gas ( |Dyson| [l9"62a|b|c| |Mehta[ |2004| ) . 

For a large class of random matrices A, the distribution 
of the eigenvalues A n can be seen as the equilibrium dis- 
tribution of fictitious point charges repealing each other 



6 



by Coulomb interaction, and submitted to an external 
one-body potential determined by the precise form of the 
probability distribution P(A). In particular, this state- 
ment is true for the Wigner-Dyson ensemble defined as 



P(A) = C 



Ne 



-l3NTrV a (A) 



(40) 



where the 'potential' V 9 (A) is arbitrary, provided the ex- 
istence of the partition function Cjy 1 . If V 9 is quadratic, 
we recover the Gaussian ensemble ([!]). To justify the 
Dyson gas picture, it is sufficient to consider the (joint) 
probability distribution of the eigenvalues (for the proof, 
see Sec. [iLC^ ): 



P({A„}) = CV-^ S({A " }) , 



N 



H 9 ({A n }) =NJ2 V 9 (^n) - E ln l A » - A ' 



(41) 
(42) 



n=l 



We recognize the Boltzmann-Gibbs distribution of a clas- 
sical gas in thermal equilibrium at a temperature T = 
1//3. The logarithmic pairwise repulsion 



V int {z) 



N 

Eh 

n=l 



A n I 



(43) 



is the Coulomb interaction in 2D, associated with the 
electric field g = (Reg, Img) represented by the resolvent 



(191: 



Ng(z = x + iy) = - V x , y V l 



N 

E 

n=l 



ReA„ y — ImA r , 



\z - A T . 



A„ 



(44) 



For Hcrmitian matrices, the Dyson gas is a two- 
dimensional Coulomb gas, experiencing the one-body 
potential V 9 , with the kinematic restriction that the 
charges move along the line ImA n = 0. 2 

The main advantage of the Dyson gas picture is that it 
allows to apply methods of statistical mechanics to calcu- 
late distributions and correlations of eigenvalues, giving 
therefore a physical intuition about the statistical prop- 
erties of eigenvalues. In particular, it is clear that the 
shape of the overall density will strongly depend on the 
one-body potential V 9 , while the correlations in the rel- 
ative positions of eigenvalues are affected by the interac- 
tion V mt and are much less sensitive to V 9 . 



Mathematically, Eq. (41) can be obtained from Eq. (WOl 



by changing variables from the matrix elements of A to 
parameters related to eigenvalues and eigenvectors of A. 
The Jacobian of the transformation contains, in particu- 
lar, a factor |V({A„})|' 3 , where 



V({A„}) = J] (A„ - A r , 



(45) 



is a Vandermonde determinant. V({A„}) is the source 
of the logarithmic repulsion in H 9 . Integrating over 
the parameters related to the eigenvectors, one obtains 
Eqs. (Ell and (|42l (iMehtal 120041). 



Equation ( 41 ) can also be proved elegantly using physi- 



cal arguments. Interpreting Eq. (40 ) as the stationary so- 
lution of a Fokker-Planck equation ( |Dyson 1972| |Mehta| 
2004 ) , it is easy to infer the associated Langevin equation 



that controls the fictitious dynamics, parametrized by 
the fictitious time r, of the independent matrix elements 
A v (t), 3 as well as the drift and diffusion coefficients of 
the matrix elements A„: 



Mi(A n ) 
M 2 (A V ) 



hm 

At->o Ar 



-N 



dV 9 {A ri 
dA n 



<AA 2 ) 



= lim 

At^o 2At 



2/3 



[1 + 6, 



(46) 
(47) 



where (. . .) denotes the ensemble average over the fic- 
titious Markov processes. This averaging must not be 
confused with averaging over different realizations of the 
random matrix A. The key point now is that we can cal- 
culate, by a second-order perturbative expansion in time, 
how the eigenvalues are modified during a time interval 
Ar : 



AA, 



AA„ 



AA 

171^71 f-l—iJ 



(p)2 
ran 

A„ 



(48) 



Averaging this equation using Eqs. (46) and (47), and 
keeping only the terms ~ O(Ar), we find (AA n ) and 
(AA 2 ), and the related drift and diffusion coefficients 
for the eigenvalues: 



Mi(A n ) = -N 



dV 9 (A n ) 
dA n 



(49) 



(50) 



2. Brownian motion of eigenvalues 

Before exploiting further the Dyson gas picture, let 
us justify Eq. (41) for P({A„}) in two different ways. 



We recognize in the drift coefficient ( 49 ) the deterministic 
force driving the point charge located at A„. In particu- 
lar, we understand in a new way the origin of the electro- 
static repulsion ( 43 ) , since in the present context it arises 



2 This kinematic res triction is suppressed for non-Hermitian ma- 
trices (see Sec. Ill I. 



•t] labels independent elements of A. Alter nati vely, we can write 



Ami, with fi ■ 



0, 



1, see Eq. 148 
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from the second-order perturbative term in Eq. (|48|) . Fi 



nally, from the coefficients ( 49 1 and ( 50 1 , it is straightfor- 



ward to reconstruct the Fokkcr-Planck equation obeyed 
by the fictitious time-dependent joint probability density 
of the eigenvalues, and its stationnary solution is pre- 
cisely the desired result (41 ). 



3. Mean-field approximation 

Once the probability distribution P({A„}) is known, 
the density of eigenvalues p(A) can formally be recovered 
by integrating it (N — 1) times. Luckily, we can avoid 
this cumbersome calculation by taking advantage of the 
Dyson gas picture. In a naive mean-field approach, the 
distribution of charges at equilibrium is found by mini- 
mizing the energy H 9 given by Eq. ([42]) . This leads to 



- d A V int {A) = Nd A V 3 (A). (51) 
Furthermore, since for Hermitian matrices ImA n = 0, 



Eqs. (43) and (44) yield 



JY 



iVRe 5 (A + ie) = ^ 



1 



A - A, 



-d A V int (A), (52) 



so that the combination of Eqs. (51) and (52) allows to 



relate Reg (A' + ie) with the one-body potential V 9 . In- 
serting the result into Eq. ( 23 ) , we obtain 



p(A) = 



1 



7rV(A-a)(fc-A) 

7T - P 



A' - A 



Let us now justify this mean-field result in a different 
way. In the large N limit, we can coarse-grain the energy 
functional H 9 : 

/oo 
AA P {A)V 9 {A) 
-oo 
j2 c roo 



Ar 
~2~ 



dAdA'p(A)p(A')ln|A-A'|. (54) 



Rigorously, when changing the integration variables from 
{A„} to the density 'field' p in the partition function, a 
Jacobian appears, which physically takes into account the 
entropy associated with p. We negl ect all corresp onding 
sub-leading terms of order In N/N (Dyson 1972). 4 The 



equilibrium of the Dyson gas is given by the extremum 
of this functional. Note that we also have to take into 
account the normalization constraint on p, which can be 



done by introducing a Lagrange multiplier c. We thus 
find: 

/>oo 

V 9 {A)- dA / p(A')ln|A-A'| + c = 0. (55) 



Differentiating Eq. (55) with respect to A we get 



dA'!^l = a A ^(A), 



(56) 



which admits the solution ( 53 ) for p defined on the com- 



pact support [a,b], as expected. 

The mean-field approach used to infer the eigenvalue 
distribution p(A) from the joint probability distribution 
P({A„}) is general and can be applied to any ensemble, 
provided that P({A„}) is known. Actually, P({A n }) can 
be found for a larger class of matrices than the Wigner- 
Dyson ensemble ( |40| . It is straightforward for any dis- 
tribution P(A) that is simply expressed in terms of the 
eigenvalues of A, e.g. through TrA or detA: P({A„}) is 
then obtained by multiplying P(A) by the Vandermonde- 
type Jacobian |V({A n })|° responsible for the logarithmic 
repulsion between eigenvalues. 



4. Application to specific random matrix ensembles 

a. Gaussian ensemble. We start by considering the Gaus- 
sian ensemble ([l]) corresponding to V 9 (x) = x 2 /A in Eq. 
(40). Using Eq. ( |53j ) with a = — b found by the normal- 
ization condition J_ b dAp(A) = 1, we readily obtain the 



(53) 

celebrated Wigner semicircle law dWigner 1955b: 



p(A) = i(4-A 

Z7T 



2\l/2 



(57) 



It states that for large N and on average, the N eigenval- 
ues lie within a finite interval [—2, 2], sometimes referred 
to as the 'Wigner sea'. Within this sea, the eigenvalue 
distribution has a semicircular form. 



b. Wishart ensemble. Our second example is the Wishart 
ensemble defined by Eqs. Q and ([5|. Let us focus on 
P(A) given by Eq. that corresponds to c = N/M < 1. 
As explained above, P({A„}) follows by multiplying Eq. 
(frl) by the Jacobian |V({A„})| 2 : 



P({A n }) = ^, M e- 2 ^« A »», 
1 N 

H 9 ({A n }) = -J2 l NA n - (M - N) In A„] 



n=l 



- ln|A„ - A m |. 



(58) 



(59) 



4 This is justified when the confining potential V s is 'strong'. |Tierz| 
l |2003| l gives explicit examples of failure of Eq. \54\ for a 'weak' 
confining potential. 



Note that if the quadratic V s is multiplied by an arbitrary con- 
stant a, the eigenvalue density is found by a simple rescaling of 
variables: Pa(A) = y/ap a —i(y/aA). 



This result has the same form as Eqs. (41 ) and (42 ), with 
the one-body potential 



V 9 (x) 







x — 









(60) 



which is repulsive in the limit x — > + . The linear and 
logarithmic contributions come from TrA and detA in 
Eq. Q , respectively. Note the difference with H entering 
in the definition of A = HH^ , for which V 9 is harmonic 
due to the term TtHH^ in Eq. ([s]) . Inserting the potential 
(l60l into Eq. (1531), we find 



p(A) = -L^-AXA-A-), (61) 
which is defined on the compact support [A_ , A+] with 

2 



A 



± 



f 



± 1 



(62) 



where the size M of the matrix T can be arbitrary, and 
in fact it will be infinite for the majority of functions 
f(ri,Tj). In Eq. ([65]), we assume N < M and C NjM (T) 



is a normalization coefficient that depends on the matrix 
T. For T = Im, we recover the Wishart case 

0- This 

shows that the eigenvalue density of the ERM associated 
with the simplest matrix T yields already a non-trivial 
result, the Marchenko-Pastur law (64 1. An explicit ex- 
ample of ERM that obeys this law for a certain range of 
parameters is given in Sec. II.I.3| For arbitrary T, in- 



ferring P({A n }) from Eq. (|65[) is a priori not easy, inas- 
much as the argument Ti(H T -1 ii^) cannot be expressed 
in terms of the eigenvalues of A. Therefore, integration 
over the independent parameters related to the eigen- 
vectors of A may be complicated. We believe, however, 
that the Dyson gas picture may be promising for ERMs, 
in particular when one is interested in more complicated 
quantities than just the density of eigenvalues, such as, 
e.g., eigenvalue correlations. 



This result was derived for c < I. It is easy to find 
the solution for c > I, by noting that, according to its 
definition (19), g is the average of 



Tr 



(JV)" 



1 rp 1 

-Tr (M) 



N-M 



z 



(63) 



where we used the cyclic permutation of the trace oper- 
ator. From Eq. (21), it is thus clear that the case c > I 
is obtained by adding N — M zero eigenvalues to p(A) . 
For arbitrary c, p(A) has the generic form 6 



p(A) 



5(A) + — v/(A+-A)+(A-A_)+, 



2ttA 



where x + = max(x,0). 



Tulino and Verdii 2004) 



The result ( 


64 


1 is the famous 


Marchenko anc 


. Pastur 


1967: 



c. Euclidean random matrices? It would be tempting to 
apply the Dyson gas picture to ERMs. This requires to 
find P({A n }) in a form similar to Eqs. ( 58 1 and ( 59 ). The 
problem is that, in order to derive P(J~A n }) with stan- 
dard tools of RMT, we need P(A), which is unfortunately 
unknow n for ERMs. However, as we discussed in Sec. 
IIA.3| an ERM A can be rewritten as A = HTH\ with 
entries Hi a that approximately behave as i.i.d. Gaus- 
sian random variables. The probability distribution of H 
is then given by Eq. (|5|). Hence, following the original 



Wishart's idea (Wishart 
the form 



irt, 



1928), we expect P(A) to be of 



P{A) = C NM (T)detA 



M—N—NTr(HT~H') 



(65) 



If Eq. jsjl is modified into P a (H) - 
ing of variables shows that p a (A) 



C N M e- aNTrHH \ arescal- 



ap a =i 



(aA). 



D. Euclidean random matrices: heuristic approach 

Before introducing advanced techniques to describe the 
spectal properties of ERMs, it might be useful to present 
simple arguments that apply for ERMs of the form Aij = 
/(r<— ij-), in the situation where the function /(r) decays 
fast enough for large r so that, in the limit V — > oo, the 
eigenvalue density p(A) depends on the density p — N/V 
only. We note that this is far from being the case for all 
ERMs; a couple of explicit counterexamples is given in 
Sec. JTLSl and ITLTl 



(64) 1. High density expansion 



For any ERM A, we can always formally write 
£f =1 ^^(k) = Ai(k)$i(k) with $i(k) = e~ ikr - and 



JV 



A i (k)=^e ik '( r 

3=1 



'/(r.-r,). 



(66) 



Suppose now that the density is large enough for the 
phase ik • (r^ — r^) to vary only weakly between neigh- 
boring points Ti and Yj. In d dimensions, this is roughly 
satisfied for p x l d ^> k. The sum in Eq. (66) can then 
be approximated by an integral, so that Aj(k) does not 
depend anymore on i, becoming an eigenvalue of A, 
A(k) = p/o(k), where 



/o(k) = J d !i r/(r) e ik 



(67) 



is the Fourier transform of f(r). This eigenvalue is asso- 



ciated with the eigenvector (e 



-ik-ri 



. . , e" 



-ik-r/, 



Sum- 



ming over the different eigenvalues labelled by k, we ob- 
tain the eigenvalue density (Mezard et al. 1999): 



P (A) 



(27T) 



^[A-p/ (k)]. 



(68) 
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For ERMs of the form Q with u = 1, Eq.([66]) is re- 
placed by 



N 

A l (k)=^[e ik '( r -- r ^-l 

3=1 

and the eigenvalue density reads 
1 f d d k 



/(r, - r,) (69) 



p(A) = :/^f[A-/, [/ (k) - / (0)]] . (70) 



2. Low density expansion 

In the low density limit p — > 0, for a rapidly decaying 
function /(r^ — i\j), the matrix elements A^ are sizable 
only if the points rj and Yj are nearest neighbors. In 
this case, the matrix A can be approximated by a block 
diagonal matrix with N/2 blocks of size 2x2. Each block 
describes the coupling between two nearest neighbors i 
and j. Its eigenvalues are A = /(0) ± /(r^ — r^). Hence, 
the eigenvalue density of A is 



p(A) 



1 



d d Ar Pnn (Ar){S[A~f(0)~f(Ar)} 



+ <5[A-/(0) + /(Ar)]}. 



(71) 



where p nn (Ar) is the probability to find two nearest 
neighbors separated by a distance Ar: 



Pnn (Ar) = dV p(Ar) d - 1 e- Vp ^ d 



(72) 



In this expression, V = Tx d l 2 /T{d/2 + 1) is the volume of 
a d-dimensional unit sphere. If f(r) is a monotonically 
decreasing function decaying to in the limit r — > oo, the 
cumulative C(A) = dA'p(A') takes a simple form: 

C(A) = ~sgn[/(0) - X}e- v Pif-^no)^\]} d + 1 (?3) 

for A <= [0,2/(0)], C(A) = 1 for A < and C(A) = 
for A > 2/(0). p(A) can be obtained as the negative 
derivative of C(A). 

For an ERM defined by Eq. ^ with u = 1, the two 
eigenvalues of the 2x2 blocks become — 2/(r, — r\,) and 
0. Therefore, the eigenvalue density reads 

p(A) = l [ d d Arp„„(Ar)<5[A + 2/(Ar)] + ^(A), (74) 



and for a monotonically decreasing function f(r) decay- 
ing to 0, the cumulative reduces to 



C(A) 



1 



-Vp{/- 1 [-A/2]} d 



(75) 



for A € [-2/(0), 0], (7(A) = 1 for A < -2/(0) and 
C(A) = for A > 0. Equation (FtS shows that half 



of the eigenvalues are zero. This unphysical result may 



come from the fact that the condition A, 



E 



cannot be applied for every block independently. A sim- 
ple way to get rid of those zero eigenvalues is to suppress 
the delta contribution as well as the prefactor 1 /2 in Eq. 
(74). The fact that half of the N points contribute to the 



remaining term can be taken into account by replacing p 
by p/2 in p nn (Ar). The corresponding cumulative 



C(A) 



-A/2]} d /2 



(76) 



was used by Amir et al. ( 2010[ ) to study the eigenvalue 
density of the exponential ERM. We review this study in 
more detail in Sec. ITlO 



E. Field representation 

In this section we discuss a field-theoretical represen- 
tation of the resolvent g{z). The starting point is Eq. 
(25), that we rewrite as 



g(z) = ~d,(hx dct(z - Ay 1 ' 2 ) 



(77) 



The determinant det(z — A) 1 / 2 can be represented as a 
canonical partition function 



Z{z) = det(z - A)- 1 ' 2 



exp 



--$ T (zI N -A)9 



, (78) 



where & T is the transpose of the vector $ = (<pi, 0jv)- 
In this way, we recast the calculation of the resolvent 
g(z) into a statistical mechanics problem of N interacting 
particles <pi with a Hamiltonian 



N N 



2 --»jvtvj< 

i=l i#i=l 

The corresponding Boltzmann-Gibbs distribution is 

1 



(79) 



P(<M) 



Z(z) 



(80) 



so that the resolvent ( 77 ) is proportional to the derivative 



of the average thermodynamic free energy, — In Z(z): 

N 



9{z) 



-^(lnZ ( z)) = 



(81) 



where 



denotes the field average with respect to 
led b 
apply 

use of the identity 



P(<$>,z) defined by Eq. (80). In order to compute 
(lnZ(z)), we apply the replica method based on a smart 



x n — 1 
In x = lim . 

n->0 n 



(82) 



The idea is to compute the right-hand side for finite and 
integer n and then perform the analytic continuation to 
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n — > 0. 7 Eq. (81) becomes 



lim-(Z"(z)) 

n-yO n 



(83) 



The quantity that we now want to evaluate is (Z n {z)) . 
It contains n copies (replicas) of the original system (|78| : 



(d^...d^)...(d^...d<^) 



exp 



2 ^ 

Q = l 



$ aJ (zljv - A) <f> a 



(84) 



The averaging (. . .) in this equation can be performed 
in different ways, depending on what we know about A. 
In the standard RMT, the averaging is performed over 
the distribution P(A) that is known. Without entering 
into details, let us describe the way to proceed with Eq. 
(84) in this case. First, we perform two algebraic manip- 



ulations: we integrate over the matrix elements (which 
is possible, in practice, for Gaussian-like distributions), 
and we introduce auxiliary fields such that integration 
over replica variables can be carried out. We thus get a 
new integral form that depends only on these new fields. 
Second, in the N — > oo limit, we find the relevant values 
of these fields by making a saddle point approximation. 
This method was originally applied to the Gaussian en- 
semble ([I]) by |Edwards and Jones (1976) who rederived 



y, it was also applied with 



the semicircle law ( 57 ). More recent! 
to Wishart matrices Q (with arbitrary covariance ma- 
trix), and the Marchenko-Pastur law (64) was recovered 
( |Sengupta and Mitra[ [l999[ ). 

For Hermitian ERMs of the form /(r^r^) = /(r, — 
r ? ), the field-theoret ical approach was first proposed by 
Mezard et al. (19991. Let us review some details of their 



calculation. For Aij = /(rj — r^), Eq. (84) becomes 



(Z n (z))cx /(d^...d^)...(d0]v...d^) 



d d ri 



d d r N 



V V 

n N 
a—1 i,j=l 



exp 



N 



a=l i=l 



* ^3 



(85) 



As explained just above, in order to perform the Gaussian 
integration over the replica fields, we introduce two sets 
of auxiliary (bosonic) fields ip a and i/j a , i.e. we insert 
into Eq. (|85| the relation 



II D hH s f 



N 



^ Q (r)-^^(r-r0 



i=l 



(86) 



7 For some models, the analytic continuation may not be unique, 
and the replica trick may fail. In a more rigorous treatment, we 
have to use the supersymmetric approach | |Efetov[ |1997[ |Haake[ 

20T0I. 



where Sf stands for the functional Dirac delta: 



6f[iP] = / D^] exp 



i / d d r V(r)V(r) 



(87) 



We then integrate out the <j> variables to obtain a field 
representation of the partition function 



(Z n (z)) = 



1 



yNn/2 



D[tp a ,^ a ]A N e Sl 



where 



r 1 " 

A= /d d rcxp - J]^(r) 2 

J L Z a=l 

n - 

5o = i E / d d r^(r)^ a (r) 

Q = l 

i n r 

+ 2 E J dr d"rV Q (r)/(r - r')V a (r')- 



(89) 



Finally, integrating out the ip fields, we get an expression 
which is a good starting point for different approxima- 
tions: 



(Z n (z)) = / D\r\e s \ 



(90) 



Si - Mn|z-"/ 2 yd d rexp _JL^^ a ( r ) 2 | 

1 - r 

+ o E / d rf rd rf rV a (r)/- x (r - r')^(r'). (91) 

Q = l 

Here / _1 is the operator inverse of / considered as an 
integral operator: 

J d d r"f-\r - r")/(r" - r') = S(r - r'). (92) 

Suppose now that we can expand the logarithmic term 
in Eq. (91). We omit terms independent of ip and apply 



the Wick rotation ip — > iip, so that the action Si becomes: 
S 1 ^pz- n ' 2 /d d rexp ^E^W 

J L Z a=l 

1 n r 

- 2 E J d'rdV^WZ-^r - r> Q (r')- (93) 



In the high density limit p = N/V — > oo, Mezard 



et al. ( 1999 ) proposed to expand the exponential term 



of the action (93), at z/p fixed. Inserting the result into 

(94) 



Eq. (90), we obtain: 



i r d a k 



pj (2n) d z - pf (k) : 
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where /o(k) is the Fourier transform of /(r), defined in 
Eq. (67 1. The corresponding density of eigenvalues ( 21 ) 



coincides with the result ( 68 ) obtained from heuristic ar- 
guments. 

In order to obtain an expression for the resolvent g(z) 
valid beyond the high density regime 



Mezard et al. 



( 1999 ) looked for the best quadratic action S v that ap- 



proximates the full problem (93): 



S v = -\\ d d rdV^ :r (r)Jif- 1 (r ) r')*(r'), (95) 

where 4 ,T = -0"). The nx n matrix K _1 (r, r') is 

obtained by minimizing the variational free energy F v = 
{S\) v — \nZ v , where Z v = / D^Je 51 ' and (—) v is defined 
with respect to the measure P v — e Slt / Z v . This leads to 
the following self-consistent equations for the resolvent 

S«: 8 



9{z) 
a(z) 



1 



a{zY 



d d k 



/o(k) 



(2ir)dl-pf (k)g(z)- 



(96) 
(97) 



These equation implicitly assume that the function /(r) 
decays fast enough for large r, so that the translational 
invariance is restored in the limit V — > oo at fixed den- 
sity p — N/V. Consequently, the resolvent g(z) and the 
density of eigenvalues p(A) depend only on the density p. 



F. Diagrammatic approach 

1. From Gaussian and Wishart ensembles to Euclidean random 
matrices 

Before discussing in details the diagrammatic treat- 
ment of Hermitian ERMs, we briefly review the results 
for Gaussian and Wishart matrices. The starting point 
of a diagrammatic computation of the resolvent ( 19 ) is 
its series expansion (27). For Gaussian-like ensembles, 



the result of averaging can be expressed through pair- 
wise contractions, such as Eq. The different terms 
(diagrams) arising from this calculation are conveniently 
collected in the self-energy <j(z) defined by Eq. (28). By 



construction, a(z) is the sum of all irreducible diagrams 
contained in the expansion of g(z), i.e. those that can- 
not be separated into two independent diagrams linked 
by the propagator 1/z. We do not detail the diagram- 
matic representation of a(z) for Gaussian and Wishart 
ensembles because these ensembles can be considered as 
special cases of ERMs, for which a diagrammatic calcu- 
lation is given below. 



8 Although Eqs. \96\ and | |97| do not appear explicitly in jMe zard 
et al. , 1999 ), it is straightforward to obtain them from the results 
presented in that paper. 



It is easy to show, using the pairwise contractions ^ 
for GUE and ^ for GOE, that the s elf-energy a{z) of 
the G aussian ensemble ([I]) is given by ( Jurkiewicz et al] 
2008I) 9 



a(z) = g{z). 



(98) 



Inserting Eq. (|98|) into Eq. (|28|), we find the resolvent 

(99) 



9{z) 



- V ' z 2 - 4 



which leads, via Eq. (21), to the semicircle law (57) 



The self-energy o(z) of Wishart random matrix ensem- 
ble Q is obtained in a similar way. The main difference 
with the Gaussian case is that we now have to distin- 
guish, when manipulating pairwise contractions ([6]), in- 
dices i = l,...,N and a = 1, . . . , M. For c = N/M < 1, 



the self-energy is (Jurkiewicz et al. 2008): 



a(z) 



1 



cl - g(z) 



(100) 



Equations «28\ and (100) lead to a quadratic equation 



for g(z), that has a normalizable solution 



2 + 1 



z-A+)(*-A_) 



(101) 



with A± given by Eq. (62). From Eq. (21), we recover 



the Marchenko-Pastur function (61) 



Historically, neither the Wigner semicircle law (57) nor 
the Marchenko-Pastur law ( |64| were derived by calcu- 
lating diagrammatically the self-energy a(z). Wigner's 



original proof (Wigner 1958) is based on an explicit cal- 
culation of the moments (A ra ) that appear in the series 
expansion (27) of the resolvent. This is somewhat sur- 



prising inasmuch as the counting procedure required to 
evaluate the moments is more complicated than the di- 
rect evaluation of the self-energy ( 98 1 . Odd moments of 



the symmetric semicircle law are zero, and even moments 
are the Catalan numbers: 



(A 2p ) = 



(2p)\ 



(102) 



P \(p + iy: 

A calculation of the Marchenko-Pastur law from its mo- 

The procedure 



ments can also be performed (Bai 1999 



is quite tricky, as we can imagine by loo 
for the moments: 



dng at the result 



(A" 



1 



n!(n- 1)! 



nh-X 



^(n + l-fc)!(n-fc)![(fc-l)!] 



(103) 



9 1 Jurkiewicz et al.\ | |2008| ) treat GUE. The GOE case is slightly 
more involved since Eq. |3jl generates two types of diagrams 
rather than one for Eq. 121. However, in the large limit, 
the second term of (J3J does not contribute to u(z) because it 
gives rise to non-p lanar d iagrams only (for the definition of these 
diagrams, see Sec. |II.F.3[ l. 
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Undoubtedly, if we are interested in the full distri- 
bution p(A), the counting procedure for evaluating the 
moments is less appropriate than the diagrammatic self- 
consistent calculation of the self-energy. The same re- 
mark holds for ERMs, as we will see shortly. 

The main object in the diagrammatic analysis of a Her- 
mitian ERM A is the operator 



Assuming translational invariance, we can eliminate one 



integral in Eq. (110), thus showing that (A™)^ 



(n) 



When two indices are equal in Eq. ( 109 ), we have a miss- 
ing factor from the sum and a missing 1/V factor from 

the average, leading to (A")<f 1} - p n ~ x . Therefore, 
(A")k has the following density expansion: 



IN N 

0(*) = (££ 
\*=i j=i 



i 

z — A 



l r i)( r i 



(104) 



For later convenience, we also define 



5k (z) = -(k|0(z)|k), 
P 



N N 



V*=l 3=1 



1 




z- A_ 


i,3 / 



(105) 



Since in the limit k — ¥ oo only terms i = j contribute 
significantly in Eq. ( 105 1, gk(z) is related to the resolvent 
PI by 



g{z) = Urn g*{z). 



(106) 



Similarly to g(z), gk(z) admits a series expansion in 
its holomorphic domain: 



, , f (A")k 



n=0 



N N 



( A ^4 EE eik ' (ri_r,,) ^ 



(107) 
(108) 



= 1 3 = 1 



2. Euclidean random matrices: high density expansion 



In this section, inspired by the work of |Grigera et al.\ 
(2011 2001b), we will present a perturbative calculation 



of the resolvent ( 105 1 by direct evaluation of moments 
{108} for ERMs of the form A ZJ = f(r h rj) = /(r< - rj). 
The moments 



N 



N 



i„+i = l 



(A") 



^(A™)« with <A»)« 



i=l 



(111) 



Let us compute explicitly the two first leading terms 

\(") _„j /\n\(n-l) 
'k 



and (A" 



in the high density regime, (A" 

We replace all functions /(r, — r 3 ) in Eq. (jl lOh by 

/ d d k/ (k)e- ik ( ri - r ^/(27i-)' i , and perform the n spatial 
integrations. For points in a box of side L and volume 

V = L d , (A n ) k n) becomes 



<A">£° = / d-kx-.d^sinc 



~ (27r)" d 
k n _i k 
2/L 



k-k, 

2/L 



2/L 



/o(k 1 ).../ (k„), 
(112) 



where, for d = 3, sine [k] = sine [k x ] sine [/c 2 '] sine [k z ]. As- 
suming that /o(k) is centered around, say, k Q , with a 
width Ak a such that k a L ^> 1 and Afc Q L <C 1, we use 
sinc[(k l - k^L/2] ~ (27r) d $(ki - kj)/L d and obtain 



<A»)W = ^/ (k)r. 



(113) 



Inserting this into Eq. (107), we obtain the crudest ap- 
proximation g^{z) for the resolvent gu.(z): 



i 



z - p/o(k) 



(114) 



This result is the 'bare' propagator that does not cap- 
ture any fluctuations of A. Using gk(z) — g^(z) leads to 
a trivial result g(z) = lim^oo g^(z) = 1/z. We thus cal- 
culate the next-order contribution (A n )| c ™ _1 ' ) , which con- 
tains two equal indices. There are two differences with 
the calculation of (A™)^ . First, we can choose the two 
positions of the equal indices. Second, for given positions 
such that we have /3+2 functions / between the two equal 
indices, 10 we replace j3 + l sine terms by (5-functions. The 
result reads 



/ ,1 1 | I 



(109) 



can be expressed as sums of n terms characterized by the 
number of repeating indices. The term with all indices 
different is 



(A 



n\{n-l) 
k 



£ b/ (k)r 



(2ir) d 



[p/o(q)] 



(/3+2) 



[p/o(k)] 7 . (115) 



{A n } (n) 



N r 



A d n d d r r , 



v V V 
x /(ri - r 2 ) . . . /(r„ - r n+1 )/(r„ + i - r x ). (110) 



10 /3 G [0, n — 2] is an integer that should not be confused with the 
symmetry index denned in Sec. |II.A.1| 



13 



Summing over n to get the corresponding resolvent ( 107 1, 
<7k(#) — 9k (z) , suppresses the restriction imposed on a, 
P and 7: 



1 



2 - p/o(k) 
1 

^ - p/o(k) 



d d q 



1 



(2ir) d z - pf (ci) 



[p/o(q)f 



(116) 



which is of the form g^(z)a (z)g^(z). The first irre- 
ducible diagram contained in the self-energy ak(z) — 
V5k( z ) — l/9k( z ) is therefore independent of k and reads 

^W = J/^^°^ a flSW- (117) 

If now we restrict the den sity expansion of the self- 
energ y to the first order (1171, <7k(z) — cr 1 (^), the resol- 
vent ( 105 ) is given by 



Sk(z) 



1 



z - p/o(k) - cr 1 ^)' 



(118) 



and, from Eqs. (21) and (106), the density of eigenvalues 
takes the form 



p(A) 



Imcr 1 (A + ie) 



[A - rW(A + ie)] 2 + paw 1 (A + ie)] 2 ' 



(119) 



For |A| > |Reo- x (A+ie)|, \lrna 1 (A+ie}\, p(A) ~ Imcr 1 (A+ 
ie)/A 2 . Using the explicit form (117) of 01, we recover 
the result (68 1. 11 This indicates that the more diagrams 



we take into account in a^z), the more accurate is p(A) 
at small |A|. Applying essentially the same procedure as 
for the calculation of cr 1 (z), it is also possible to compute 
higher-order contributions a % v (i > 1) of order l/p l to 
the self-energy a^(z) = cr 1 (z) + a\{z) + . . . , even though 
the combinatorial rules presented in the recent literature 



( Ganter and Schirmacher| 2011| Grigera et al. 2011 ) are 
quite involved. A simple way to improve Eq. ( |117[ ) is 
to replace the bare propagator g^(z) in the integrand by 
9q(z), thus obtaining a system of self-consistent equa- 
tions for <r(z) and gu(z). 

With minor modifications, the analysis of this section 
can be repeated for ERMs defined by Eq. ^ with u — 
1. An important complication in this case is that the 
self-energy a 1 becomes k-dependent. The self-consistent 
equations for g-^{z) and a^z) read (Grigera et al. 2001a|: 



9k(z) 
0k(z) 



*-p[/o(k) 
1 f d rf q 
PJ (2tt) 



/o(0)] 
l ^{p[/o(q)- 



(120) 



/ (k-q)]} 2 ffq (z).(121) 



11 Another way to recover Eq. j68| l is to compute the series j27| 
with (A") ~ <A n )( n ) calculated following the same procedure as 
for (A™)'"'. 



(a) 

H = J 



H. T nHl 

ia ap pj 



\\ f // 



(b) 



= 1/JV 
= Nx, /Tfj = Tri" 



FIG. 1 (a) Diagrammatic representations of the matrices H, 
H\ and A = HTHK Full and dashe d lines propagate in 
the bases {r^} and {i/j a } defined in Sec. II. A. 3 respectively; 
T = pA. (b) Diagrammatic notation for pairwise contractions 
Q and loop diagrams for any scalar x in the basis {r^} and 
for any operator X in an arbitrary basis {ip a }- 




FIG. 2 Diagrammatic expansion of the resolvent g(z). A 
horizontal straight line reprepresents the propagator 1/z. 



3. Euclidean random matrices: self-consistent equations 



In this section, we derive self-consistent equations for 



the operator (104), using the representation (10) for an 

(rj|A|rj). We recall that the 



ERM Ai 



matrix H is random but independent of the function /, 
whereas the matrix T depends on / but is not random 
(see Sec.|II.A.3l. 



a. Diagrammatic calculation. We start by expanding the 
resolvent (19) in series in 1/z: 



!{*) = 



Tr 



'-A 1 - 

z z 



1 ,1 A 

-A-A- 

z z z 



(122) 



where averaging (...) is performed over the ensemble of 
matrices H . As explained in Sec. |II.A.3[ we assume that 
H has i.i.d. complex entries distributed according to the 
Gaussian law |5]). Using the properties of Gaussian ran- 
dom variables (such as the Wick's theorem) , the result of 



averaging in Eq. ( 122 ) can be expressed through pairwise 
contractions To evaluate efficiently the weight of dif- 
ferent terms that arise in the calculation, it is convenient 
to introduce diagrammatic notations. The matrices H, 
H\ and A will be represented as shown in Fig. [lja). 

Each contraction ^ brings a factor 1/N, and each 
loop corresponding to taking the trace of a matrix brings 
a factor N, see Fig. [ljb). In the limit N — >• 00, only 
the diagrams that contain as many loops as contractions 
will survive. These diagrams are those where full and 
dashed lines do not cross. Therefore, the leading order 
expansion of the resolvent ( |122 1 involves only diagrams 
which are planar and look like rainbows, see Fig. [2] where 
we show the beginning of the expansion of g(z). Note 
that the prefactor 1/N of Eq. (122) does not appear in 
Fig. [2] because it is compensated by the external trace. 
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_ Tr( T ) 

2 5 TV 3 



FIG. 3 A typical non-planar diagram appearing in the expan- 
sion of the resolvent g(z). Its value follows after application 
of 'Feynman' rules defined in Fig.[TJb). It does not survive in 
the limit TV — > oo. 




FIG. 4 Diagrammatic expansion of the self-energy a(z). 
Braces with arrows denote parts of diagrams that are the be- 
ginning of the diagrammatic expansion of the resolvent g(z). 



An example of a non-planar diagram is represented in 
Fig. [3] It vanishes in the limit TV — > oo. 

The self-energy o~(z) is the sum of all one-particle irre- 
ducible diagrams contained in zg(z)z. The first dominant 
terms that appear in the expansion of o(z) are repre- 
sented in Fig. |4j Under a pairwise contraction, we rec- 
ognize g(z) depicted in Fig. [2j After summation of all 
planar rainbow diagrams in the expansion of Fig. [4] and 
application of 'Feynman' rules defined in Fig. |Jb), the 
self-energy becomes 



*(z) = ^Tr 



TrT 



T 

1 - g(z)f_ 

N l 



(123) 
(124) 



where f = pA, and Tr denotes the trace of an operator. 
Inserting Eq. (1231 into Eq. (28), we obtain: 



1 



i 

TV" 



Tr 



T 

g(*)T 



that is a closed equation for g(z). Noting that 12 
Trf = pTrA = (Tr N A) = N(A), 



(125) 



(126) 



we conclude that TrT/TV in Eq. (124 1 leads to a shift in 
the distribution of eigenvalues p(A). 



Before commenting on the result (123), let us see how 



the operator (104) can be expressed through the solution 



g(z) and f. In the basis {ip a }, Eq. (104) reads 



'-A 1 - 

z z 



N N 
1 = 1 ] = 1 

1 A A 

- -A- A- + . . 

z z z 



H jP , (127) 



where we used the definition ( 11 1 of the matrix H. In Fig. 
[Sj we represent the beginning of the expansion of O a p / p 
with the diagrammatic notations of Fig. [ija) . Note that 
all diagrams in Fig. [5] are irreducible. As it was the case 
for <j(z), we recognize the expansion of g(z) under pair- 
wise contractions. After summation of planar diagrams, 
the operator O(z) is finally given by: 



6{z) 



pgjz) = p 

1 — g(z)T z — T — <r(z) 



(128) 



b. Low density limit. Without loss of generality, let us 
assume that the diagonal elements of the matrix A are 
all equal, An = (A). At low density p — > 0, an ap- 
proximation of the self-energy ( 123 ) can be obtained by 



neglecting the term g(z)T in the denominator: 
a(z) ~ + T ^ 2 ^g(z) = (A) + VarA. 9 (z). (129) 



N 



N 



The last equality of Eq. (129) follows from 



Tr(f 2 )=p 2 J J d d rdV|/(r,r')f 

(Tr^)) 2 



(Tr/v(A )) 



N 



NVarA. (130) 



The implication of Eq. ( 129 ) is that in the limit of low 
density p, the approach described in this section yields 
for all ERMs an eigenvalue density identical to that of a 
Gaussian matrix: 



p(A) 



1 



2?rVarA 



V^VarA- (A - (A)) 2 , (131) 



where the variance VarA is given by Eq. (130). This re- 



sult holds for ERMs for which the representation ( 10 ) 
with matrices H having i.i.d. elements and T having 
a sufficiently large number of non-zero eigenvalues, is a 
good approximation. When considering specific exam- 
ples of ERMs in Sec. we will see that both matrices 
that fall into this category (Sec. II.I.3 and 11.1.4) as well 
as those that do not (Sec. II.I.1[ ) exist. The low-density 
limit of the latter can be understood using the heuristic 
analysis of Sec. |ILD| 



12 From here on, we denote by Trjv the trace of a TV X N matrix 
when confusion is possible with the trace of an operator. 



c. Relation to previous results. Finally, let us show how 
the various approximations for g{z), o~(z), and gk(z) 
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< a .''/J'\«> + ( a ilA\.±JlA\g) + (JLA\..T._<LAi..f.JLA\s) + 




= a ,'£r5v, a 8 ap + a £g&..T../£g$\p + J^\..f..!t^\..f.J&.p + ... 

FIG. 5 Diagrammatic expansion of O a p/p- Braces with arrows denote parts of diagrams that are the beginning of the 
expansion depicted in Fig. [2] 



derived in the previous sections can be recovered from 
Eqs. (1231 and (128 1. We need to assume that 



/(k,k') = (k|A|k' 
1 

_ V 



d d r d d r'e- i(k r - k ' r,) /(r - r') (132) 



is diagonal: /(k, k') ~ (k|^i|k)(5kk' = /(k)<5kk', which 
is not exact in a finite volume V. In the momentum 
representation, Eqs. ( 123[ ) and (1051 read now 



a{z) 



d d k 



/(k) 



-(k|0(z)|k) = g k (z) ~ j± rT , 

P z- pf(k) - a{z) 



(133) 
(134) 



can be further approximated by 
Hence, 



where /(k) = (k|A|k 
/o(k) defined in Eq. 
identical to Eq. (97). If the integrand of Eq. (1 1331) is ex- 



panded in series in p, Eq. (134) becomes consistent with 
Eqs. (117) and (118). This means that the approximate 



Eq. (133) becomes 



self-energy (117) corresponds to a truncation of the ex- 
pansion depicted in Fig. [4] after the second diagram. 



d. Solving self-consistent equations in practice. The solu- 



tion of Eq. ( 125 ) for a given matrix A can be greatly 



facilitated by a suitable choice of the basis in which the 
trace appearing in this equation is expressed. In addition 
to {r} and {k a }, a basis of eigenvectors \TZ a ) of T can 
be quite convenient. The eigenvector \H a ) obeys 

(r\f\K a ) = p [ dV/(r,r')K a (r') = (i a K a {r), (135) 
Jv 

where fi a is the eigenvalue cor respo nding to the eigen- 
vector \lZ a ). In this basis, Eq. (125) becomes 



1 1 



Pa 



(136) 



G. Free probability theory 

The term 'free probability theor y' designa tes a disci- 
pline founded by Voiculescu (1983) [see also (Voiculescu 



et al. 1992 )] in order to solve the following problem: can 



we say anything about the spectral properties of the sum 
of two matrices X\ + Xi when the spectral properties of 
the summands X\ and X 2 are known? Unless the two 
matrices commute, knowing their eigenvalues is, in gen- 
eral, not enough to find the eigenvalues of the sum. How- 
ever, the free probability theory identifies a certain suffi- 
cient condition, called asymptotic freeness, under which 
this problem can be tackled without involving the eigen- 
vectors of the matrices. The notion of asymptotic free- 
ness is equivalent to the notion of statistical indepen- 
dence that we are familiar with for random variables. It 
is a generalization of the latter to the case where the 
variables — here, the matrices — do not commute. 



1. Theoretical framework 

Let us briefly recall basic properties of independent 
variables. We denote by p x the probability density of 
the variable x, by g x (z) = (e zx ) = J^n^o z "/ n! 
its characteristic function, and by r x (z) = h\g x (z) — 
J2 n >o c x.nZ n its cumulant generating function. For two 
independent real random variables x\ and X2, the follow- 
ing relations hold: 



Pxi+x 2 = Pxi * Px 2 i 



r 



Xi+X 2 



= r 



Xl 



' X 2 ) 



(137) 
(138) 
(139) 



where V denotes the convolution. We will see that these 
relations find their equivalents for asymptotically free 
matrices. 

By definition, two Hermitian matrices X\ and X^ are 
asymptotically free if for alH € N and for all polynomials 
Pi and qi (1 < i < I), we have (Tulino and Verdii 2004[ ) 



( Pl (X 1 )) A = (q l (X 2 )) A = 

=> (p 1 (X 1 )q 1 (X 2 )...p l (X 1 )q l (X 2 )) A 



(140) 
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where the expectation value 



is defined as 



then the 5-transform is 



(X), 



N 



(TrX) 



(141) 



The interpretation of the formal definition ( 140 1 is the 
following: two matrices are asymptotically free if their 
eigenbases are related to one another by a random rota- 
tion, or said differently, if their eigenvectors are almost 
surely orthogonal. 



From the definition ( 140 ), it is easy to compute various 
mixed moments of X\ and X 2 . Considering Xi = Xi — 
(Xj) A that obey (^i)a = (A 2 )a = 0, we obtain from 



Eq. (140) 



<XiX 2 ) A = (X 1 ) A (X 2 ) A 



(142) 



Note that this last condition is not enough to define 
asymptotic freeness, since matrices do not commute. For 
example, from Eq. (140), forth moments read 



(X 1 X 1 X 2 X 2 ) A = (x 2 ) A (xf > A , 
(M 2 M 2 ) A = (x 2 1 ) A (x 2 f A + (x^l (xf >, 
-{x 1 )l(x 2 ) 2 A . 



(143) 



Free cumulants are defined such that the sum prop- 
erty ( 139 1 is preserved for the generating function of the 



free cumulants, the so-called Tvl-transform (Jarosz and 



Nowak[ |2006| |Tulino and Verdu| |2004[ ). Interestingly, 



the 7£-transform is simply related to the Blue function 
(29) — the functional inverse of the resolvent g(z) — by 
Eq. (30). 13 The 7?.-transform of the sum of two asymp- 



totically free matrices X\ and X 2 obeys: 

TZ Xl +x 2 (z) = K Xl {z) +TZx 2 (z). 



(144) 



Hence, the problem of finding the eigenvalue distribution 
of the sum of two free random matrices is straightfor- 
ward. Applying successively Eqs. (29), (30), and (144), 



one readily infers gx x +x 2 from gx ± and gx 2 - The steps 
of the algorithm are as follows: 



gXi -> Bx z -> Tlx, -> K Xl +x 2 
->• B Xl +x 2 gx 1 +x 2 - 



(145) 



There is an analogous result for the product of free ma- 
trices, which involves the so-called 6> -transform (Tulino 



and Verdu 2004). If we define x( z ) as a solution of 



1 

xO) 



1 



l = z, 



(146) 



13 Note that g(z) plays the role of a free characteristic function, 
see Eqs. ^26\ and \27\ . | Jarosz and Nowak| ( |2006| ) provide an in- 
sightful discussion about the relation between the free cumulants 
and the moments (A"). 



S(z) 



1 + z 



(147) 



Equations (146) and (147) are equivalent to the following 



implicit equation for S(z): 

S{z)K[zS(z)} = 1. (148) 
The 5-transform of the product of two asymptotically 



free matrices Xi and X 2 satisfies (Tulino and Verdii 



2004) 



S 



(z) = S Xl {z)Sx 2 {z)- 



(149) 



Therefore, the S'-transform plays a role analogous to the 
7£-transform for products (instead of sums) of free ma- 
trices. The recipe to find the eigenvalue density of X\X 2 
is analogous to Eq. ( 145 ) : 



9x % -> Xx, -> Sxi -> S Xl x 2 -> XXtX 2 -> gxtX 2 - (150) 



2. Gaussian and Wishart ensembles revisited 

A good attitude when searching for the eigenvalue den- 
sity of a given matrix, is to look for a possible decompo- 
sition of the latter in a sum or product of free matrices 
for which resolvents are known. Let us apply this idea 
to recover in a new way the now familiar semicircle and 
Marchenko-Pastur laws. 

Let us first consider a matrix A from the Gaussian 
orthogonal ensemble (GOE), with the pr obab ility distri- 
bution P(A) = CV-t Tr(A 2 )^ From Eq ( ji38| ) , it is clear 
that the distribution of the variable x\+x 2 , where X\ and 
x 2 are independent Gaussian random variables with vari- 
ances cr 2 , is still Gaussian with a variance 2a 2 . We can 
therefore decompose any Gaussian matrix A in a sum of 
two independent rescaled matrices A\ and A 2 that have 



the same probability density P: A 



\/2 



{A 1 + A 2 ). In 



addition, two independent Gaussian matrices are asymp- 
totically free. Indeed, since the measure P(A) is invariant 
under orthogonal transformation, rotation matrices Oi 
and 2 , diagonalizing Ai and A 2 , respectively, are ran- 
dom rotations over the orthogonal group. This means 
that the rotation 0[0 2 from the eigenbasis of A\ to that 
of A 2 is also random, which is precisely the intuitive def- 
inition of asymptotic freeness iTulino and Verdu (2004) 



give a formal proof of this statement] . The additive prop- 
erty of the 7?.-transform and the scaling property ( 33 ) 
yield: 

TZa(z) =Ka x (z)+11a 2 {z) = V2K A ( A= J ■ (151) 
V2 Vv2/ 

A so luti on of this equation is TZ^(z) oc z. According to 
Eq. {fr}, K'(0) = (A 2 ) = (TrA 2 )/N = 1, so that 



K(z) 



(152) 
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As could be expected from Eqs. (31 1 and (98), this is 
the 72-transform of the semicircle law. Thus, we can 
claim that the semicircle law is the free counterpart of the 
Gaussian distribution in the classical probability theory. 

In order to use the powerful arsenal of free probability 
for Wishart matrices, we decompose the A x A matrix 
A = HH^ as 

M 

HH^ = J2 h (a)t h(°» with = (H* a ,...,H* Na ). 

a=l 

(153) 

The spectrum of each matrix h^ Q ^h^ Q ^ is simple because 
it has only one non-zero eigenvalue A a — ||h( Q )|| 2 — 
X!j=i \Hia\ 2 , associated with an eigenvector h' a '*, The 
(A — 1) other eigenvectors associated with zero eigen- 
value form the basis in the hyperplane perpendicular to 
the vector h^* . Since the vectors h( a ) are uncorrelated, 
we can replace the resolvent of the matrix h^^h^) by 



l 

N 



N- 1 



1 



(154) 



where we used (An) = 1 ((|_ffi Q | 2 ) = 1/A). Inverting this 
relation gives 



72 h 



1 

2z 
1 

A 1 



(z-l) 2 + — 



O 



1 

A 2 



(155) 



For independent vectors hS a \ that have independent en- 
tries with variances equal to 1/A and identical means, it 
can be shown that the matrices hwthW are asymptot- 
ically free ( |Tulino and Verdu||2004| | . Thus, 



M 

72/mt (z) = ^2 ^h(°)th(°) (2) (156) 

a=l 
1 1 



C 1 — Z 



(157) 



where c = N/M. This is the 72-transform of the 
Marchenko-Pastur law, as could also be obtained from 
Eqs. (31) and ( |100 ). It is interesting to note that if 
we took the Ath classical convolution (by inverting the 
sum of cumulant-generating functions) of the distribu- 
tions of the variables A a , we would obtain asymptotically 
(A, M — > 00 at fixed c = N/M) the Poisson distribu- 
tion. However, the distribution that we obtain by taking 
the Ath free convolution (by inverting the sum of 72- 
transforms) is the Marchenko-Pastur law. The latter is 
therefore the free analog of the Poisson law in the classi- 
cal probability theory ( Tulino and Verdu , 2004 1 . Another 



simple proof of this law was given by 
using the product decomposition of t 
the properties of the 5-transform. 



Janik et al. 



(1997a 



re matrix HH J and 



3. Application to Euclidean random matrices 

By analogy with results for Wishart ensemble, it is 
straightforward to apply the toolbox of free probabil- 
ity to ERMs. We start with the decomposition A = 
HTW , where the basis {ip a }, that defines Hi a through 
Eq. (11), is assumed to be the eigenbasis of the opera- 
tor T: T\tf> a ) = Hafyct)- The matrix A is conveniently 
rewritten as 



M 



(158) 



where h( Q ) is defined in Eq. (153). As explained above, 
the M matrices h^^h^) are asymptotically free, as long 
as the vectors h( Q ) are independent. Hence, 



M 



K HT Hi{z) = X^oh<=>th<°>M 
a=l 
M 

= ^ Ma^hl°)th(°) (Pa z ) 
a=l 

1 M 

1 >. ^ jj, a 
a=l 



Itt 

A 



A I 



T 



1 - zT 



1 




= — Tr 


A 




1 


"1 


cz 


z 



T 



1 - zT 
1 



-9t 



- 1 



(159) 

(160) 
(161) 

(162) 



Equation ( 159 1 follows from Eqs. ( 144 ) and ( 33 ), whereas 



Eq. (160) — from Eq. (155), and Eq. ( 162 1 — from gr(z) = 

J2 a =i l/( z — ^a)M- Using the definition ( 147 ) of the <S- 
transform, one also easily shows that Eq. ( |162 ) is equiv- 
alent to 



S 



HTHf 



1 



z + 1/c 



S T {cz). 



(163) 



For completeness, we now derive Eq. ( 163 ) from Eq 



( 149 1 . From the definitions of the resolvent g(z) and the 
5-transform, one can check that, for arbitrary matrices 
A and B of size A x M and M x A, respectively, 



Sab{z) 



1 



S B a{cz). 



(164) 



z + 1/c 

Applying this result to A = HT and B = H\ we obtain 
z+1 



■5/iTfft( z ) — — , / S H >, HT {cz), 
z + 1/c 



z+1 
z + 1/c 



S H \ H {cz)S T {cz). 



(165) 



Here we made use of the fact that the deterministic 
matrix T and the random matrix H are asymptoti- 



cally free. Besides, the combination of Eq. ( 148 ) with 
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fttftfffcO = cK HH i(z/c) = c/(c-z) gives 



c + z 



(166) 



From Eqs. (165 1 and (1661, we finally recover Eq. (163) 



The result (1631, or equivalently its operator form 



( 161 1, is in perfect agreement with the solution obtained 
by a diagrammatic approach in Sec. |II.F.3| Indeed, the 
self-energy a(z) = lZ[g(z)] inferred from Eq. (161) is ex- 



actly the result ( 123 ). It is worth recalling that Eq. ( 161 ) 



was obtained from the asymptotic freeness of the matri- 
ces hWhW, that holds as long as the elements H ia are 
i.i.d. with a finite second moment ( |Tulino and Verdu 
2004). In particular, it means that Eq. (123) is valid 
even if Hi a are not Gaussian variables. Therefore, Gaus- 
sian hypothesis, t hat lar gely simplified diagrammatic cal- 
culations in Sec. II. F. 3 turns out to be not essential. 14 



In particular, this remark holds for the Wigner semicir- 
cle and the Marchenko-Pastur laws, 15 and justifies their 
large degree of universality. 

As far as ERMs are concerned, we conclude that the 
only assumption that may limit the applicability of Eq. 
(123) is the independence of the vectors h( a ). We know 
that their covariance matrix is proportional to the iden- 



tity (see Sec. II. A. 3 ), but this is not enough to insure their 
independence, precisely because Hi a are not Gaussian 
random variables. We will see how correlations between 
Hi a limit the applicability of Eq. ( 123 ) by considering 
specific examples of ERMs in Sec. 



H. Renormalization group approach 

Renormalization group (RG) approaches are powerful 
methods to identify phase transitions and study univer- 
sal properties in their vicinity. Most RG procedures that 
are applied to both ordered and disordered systems renor- 
malize the space in a homogenous way. This choice, nat- 
ural for translationally invariant systems, is not the only 
possible one when w e deal with disordered systems ( Igloi 
and Monthus 2005). Indeed, for the latter, space can 



also be renormalized in an inhomogeneous way in order 
to 'adapt' to local fluctuations due to disorder. Such real- 
space RG was first developed in the context of quantum 



spin chains (Dasgupta and Ma 1980 Fisher[ 1995 Ma 



et al. 1979 1 and then was successfully applied to various 
quantum and classi cal systems; see the revie w by |Igloi 
and Monthus! (120051) for works before 2004 and (IMonthus 



and Garel 2011 ) for more recent references. 



To illustrate the application of RG approach in ERM 
theory, let us consider a matrix ensemble defined by Eq. 



Q with u = 1 and /(r, — Vj) a positive function de- 
caying to at large distances. Follow ing |Amir et al. 
(2010) and Monthus and Garel (20111, we analyze the 



eigenproblem AR n = A n R n with the RG approach mak- 
ing use of the following mechanical analogy. Imagine 
a network of unit masses connected by springs, such 
that Aij is the random spring constant between masses 
i =/= j . Then, the displacement Xi of the mass i obeys 
d 2 Xi/dt 2 — J2j AijXj. Therefore, the eigenvalues A„ of A 
are related to the frequencies uj n of vibrational modes by 
A n = — The frequencies uj n can be found by eliminat- 
ing iteratively the modes with highest frequencies. Ini- 
tially, if the density is low enough, the highest frequency 
corresponds to the vibrational mode associated with the 
pair of the two closest masses: w = [max( J 4 i j)//i] 1/ ' 2 , 
where /i = mim,2/(mi + •m^) = 1/2 is the reduced mass. 
Provided that uj is high, this mode will not be coupled 
to the lower-frequency modes of the system and can be 
decimated: the relative displacement of particles i and j 
is set to and the pair is assumed to behave as a single 
mass rriQ = mi + rri2 connected by springs to the N — 2 
other masses a, with spring constants Aq u = A\ a + A 2a . 
This decimation procedure can be repeated iteratively. 
At step k, where we have N — k + 1 masses m, connected 
by springs of strength A^ , the mode associated with the 
renormalized frequency scale 



Aij(l/m,i + I /fhj), 



(167) 



is decimated and the two corresponding masses m, and 
fhj Q are replaced by a single mass 



fh G = fh io + fhj 
with the new spring constants 



.4 



Go 



Ai 



■At 



(168) 



(169) 



This renormalization scheme is consistent only if the new 
generated frequencies at each step are smaller than the 
decimated frequency, so that ui decreases along the RG 
Sow. This restricts the analysis to the regime of low 
density in which the distribution of couplings Aij is suffi- 
ciently broad ([AmiF^F^L 2010 Fisher 1995). It is also 
worth noting that alternative or more sophisticated RG 
rules can be considered (Monthus and Garel 2011). In 



particular, such rules may be needed to take into account 
the modes associated with individual masses and follow- 
ing the slow motion of their neighbors adiabatically. Such 
modes are neglected in the RG scheme described above. 



As the RG flow progresses, the eigenvalues A = 



of 



A are obtained in an ascending order in the range from 
—2/(0) to 0. After the step k, the flow has reached the 
typical value A given by the relation 



14 A rigorous diagrammatic proof of Eq. \123\ with the finiteness of 
the second moment of Hi a as the only assumption is nontrivial. 

15 We are not aware of diagrammatic proofs of the Wigner semi- 
circle and the Marchenko-Pastur laws that would not invoke the 
Gaussian assumption. 



k = N dA'p(A'). 



(170) 



On the other hand, the typical size m c of the eigenvector 
associated with A is given by the the typical mass of the 
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clusters of points at step k. By construction, this mass 
grows along the flow, starting from pairs of points for the 
smallest eigenvalue and converging to clusters of N points 
for the largest one. If we assume that the formation of 
clusters is a random process independent of the choice of 
the frequency scale ( 167 1 , then it is easy to verify that 

m c = N/(N - k) after the step k ( |Amir et al~] |2010[ ). ( |Gray| |2006| |Grenander and Szegoj |1958 ) 



If now we discretize Eq. ( 172 1 on a grid of M equidis- 
tant points with a spacing Ax, we see that A Qi = p ai j Ax 
are the eigenvalues of the Euclidean M x M matrix 
B with elements By = f(iAx,jAx) = 
j\ 2 /2a 2 ). The matrices with elements Bij 



on the difference 



Combining this result with Eq. (170) yields 



exp(— Ax \i — 
that depend 
j only are called Toeplitz matrices 
Of the whole 



m c (A) 



1 



C(A)' 



(171) 



where C(A) = J^dzp(z) is the cumulative eigenvalue 
density. In the regime of low density, it is given by Eq. 
pi 



The RG rules ( |167[), (|168[ ) and ( |169[ ) are explicitly im- 
plemented in Sec. |ILI.2| where we apply them to study 
the eigenvalue density of an ERM defined through an ex- 
ponentially decaying function /(r^ — Tj), We stress that 
these rules only apply for ERMs defined by Eq. ( 9]) with 
u = 1 and cannot be used at u = because in this case 
the spring analogy does not hold anymore. 



I. Application of the general theory to specific Euclidean 
random matrix ensembles 

1. Gaussian Euclidean random matrix 

Historically, the Euclidean matrix defined through the 
Gaussian function /(r^r,-) = /(r-j — r\j) = exp(— — 
Yj\ 2 /2a 2 ) was the first to be considered in the litera- 



ture ( |Mezard et al.\ 1999). Assume that the N points 
that determine the particular realization of our ran- 
dom matrix are randomly chosen in a cube of side 
L a centered at the origin. In the low density limit 
pa 3 — > (with p = N/V and V = L 3 ), we can apply 
heuristic arguments of Sec. |H.D] and calculate the eigen- 
value density by differentiating Eq. ( 73 ) . This yields a 
sharp peak at A = 1. To obtain p[K) at higher den- 
sities, we will follow the approach developed in Sec. 
II.F.3.d First, we analyze Eq. ( |135[ ) for the eigenval- 
ues p a and eigenvectors TZ a (r) of the operator T. Be- 
cause the function /(ri,rj) is separable, i.e. it can be ex- 
pressed as f{ri,rj) = /(x;, x J )f(y i ,yj)f(z l , z 3 ), p a and 



7£ a (r) obeying Eq. |135|) can be written as products of 
corresponding eigenvalues p Ui and eigenvectors TZ ai (r) 
of a one-dimensional problem: p a = pp ai p aj p ak and 
TZ a (r) = TZ ai (x)7i aj (y)7l ak (z), where p a% and 72. Qi (r) 
obey 



L/2 



dx'f(x, x')TZ ai (x') = p a Jl a% (x) (172) 

L/2 

and the index a has now three components: a — 
{ai, ctj, afc}. Equation (136) for the resolvent g(z) can 
be then rewritten as 



— — V- 

'z) + N ^ 1 



(173) 



arsenal of powerful analytical tools developed for Toeplitz 
matrices, we will make use of the so-called fundamental 
eigenvalue distribution theorem of Szego: under suitable 
assumptions and for M — > oo, the average of F(X a .) over 
all eigenvalues \ ai converges to the average of 
over £ £ [0, 2w), where 



/(0= E /(Aafc)exp(ifcO. 



(174) 



k=-c 



Expressing p ai through X ai in Eq. ( 173 1 and applying the 



above theorem to the sum over eigenvalues, we obtain 



1 



g(z) (2ir) 3 N 

'• Jff dfrdfrdfr p/MlMlM 

o i-AxW(a)/(6)/te)sW 



(175) 



Without 
now that 
ing to an 
and then 
the limit 
ous limit 
equation 



giving technical details of derivations, we note 
the series in can be summed up, lead- 

expression for /(£) involving special functions, 
the integrals in Eq. (175) can be carried out in 
of Ax — > 0, which corresponds to the continu- 
that we had initially in Eq. ( 172 ). The resulting 
for the resolvent g{z) 



is 



1 



g(z) (27rf/ 2 pa 3 g(z) 



U 3/2 [(2ir)V 2 pa 3 g(z)}, (176) 



where Li 3 / 2 ( z) is the polylogarithm function (Prudnikov 



et al. 19901. An important feature of this equation is 



that it depends on a single parameter characterizing the 
system — the dimensionless density of points pa 3 — but is 
independent of the total number of points N and the 
system size L. This is a consequence of the fast decay of 
f(ri,Tj) with |rj — Tj\ in conjunction with the assumption 
of L > a. 

In the limit of low density pa 3 — > 0, the solution of Eq. 
( fl76| is g(z) = (z- l)" 1 and hence p(A) = 5JA - 1). A 
better approximation is obtained from Eq. (73). When 
the density pa 3 is increased, the eigenvalue distribution 
widens, but a narrow peak at A = 1 resists much longer 



than predicted by Eq. ( 176 1, as can be seen from Fig. |6j 
where we compare the eigenvalue density following from 
Eq. ( 176 ) with the results of numerical diagonalization at 
pa = 0.1 and three different numbers of points N. Even 
though the height of the peak at A = 1 decreases with 
N, it is clear that the theor y fails to describe p(A) in this 
regime (Mezard et al. . 1999 1. At higher densities pa 3 > 1, 



the peak disappears and a satisfactory approximation to 
p(A) is obtained by taking the high-density limit of Eq. 
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FIG. 6 Eigenvalue density of TV x TV Euclidean matrix A 
with elements Ay = exp(— \n — | 2 /2a 2 ) at an intermediate 
density of points pa? = 0.1. The analytic result obtained 
by solving Eq. ( |176[ ) (dashed line) is compared to numerical 
diagonalization (solid lines) . The inset is a zoom onto a part 
of the main plot around A = 1. 
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FIG. 7 Eigenvalue density of TV x TV Euclidean matrix A 
with elements Aij — cxp(— ||rj — r_j ] 2 ) at a high density of 



points pa? = 1. The approximate Eq. (1771 (dashed line) is 
compared to numerical diagonalization (solid lines). 



( 176 1 that is very close to the result obtained by Mezard 
et al.\ (Jl999 ) in the high-density approximation [see also 



Sec. |II.D.l| and, in particular, Eq. (I68]) 1 ■ 1 5 ' 

(2ttW 2 



p(A) 



1 



V27r 2 pa 3 A 



In 



pa 



,1/2 



A 



(177) 



We compare this equation with the results of numerical 
simulations in Fig. [7| 



In conclusion, the heuristic analysis of Sec. II. D cor- 
rectly describes p(A) in the limit of low density, when it 
is peaked around A = 1. In contrast, the self-consistent 
equations of Sec. |II.F.3| apply at higher densities. None 
of the two approaches correctly follow the transition be- 
tween the two regimes. In particular, the reason of failure 
of Eq. ( 176 ) resides in the assumption of independence 
of elements Hi a of the matrix H in the representation 
( 10 1 of an arbitrary ERM. We checked that if, instead 
of taking Hi a = exp(zk Q rj)/-\/rV, we generate random 
matrices A = HTH^ numerically using matrices H with 
independent elements, 17 the eigenvalue density of result- 
ing random matrices nicely follows Eq. ( 176 1. As we will 



see in the following, the assumption of mutual indepen- 
dence of Hi a will limit the applicability of the analysis 
relying on Eq. (1251 to other ERMs too, although in a 



less important way. 



2. Exponential Euclidean random matrix 

Another example of ERM defined using a rapidly de- 
caying function /(r,; — r^) can be obtained by taking 
ffa—Tj) = exp(— |r<— Tj\/£). Let us study the eigenvalue 
density of the matrix A defined by this function and Eq. 
([9]) with u = 1. Such matrices appear in various contexts; 
as an example, in Sec. |iV.B"1 we will use them to study the 
electron glass dynamics. If we assume that A^ represent 
the hopping rates between randomly distributed, expo- 
nentially localized states i, the probability pi to occupy 
the state i obeys the master equation Api /At = ^2 ^ijPj ■ 
We emphasize that the symmetry constraint ■ A^ = 
combined with the fact that A^ > has a profound im- 
pact on the eigenvalues A of A: they are all negative, 
whereas they would be positive is the diagonal elements 

= 1 [i.e., for u — in 



Eq. im] (Mezard et al. 



A 


i — / 




1999) 



Amir et al. (2010) studied the eigenvalue density p(A) 
of the exponential ERM A in the limit of low density 
p£ d <C 1. In the space of arbitrary dimensionality d, their 
main result can be rederived using the heuristic approach 
of Sec. II. D. 2 Inside its support [—2,0], p(A) is readily 



found by differentiating Eq. ( 76 ) 



, . s -dVpi d \-\n{-A/2)] d - 1 e- v P^- x <- K l 2 ^ d / 2 
KA) = 2A 

(178) 

This equation shows that p(A) ~ — 1/A for all d over 
a broad range of A's. We compare it with the result 
of numerical diagonalization in Fig. [8] for d — 3. Good 
agreement is found as long as p£ 3 < 10 -3 . At larger den- 
sities, a more sophisticated theory is required. This the- 
ory should include diagrams similar to those taken into 



17 We tried different statistical distributions of Hi a , like, e.g., the 
16 This approximation does not describe the rounding of p(A) at circular Gaussian distribution and Hi a = exp(itpi a )/V~N , with 

small A and does not obey the normalization of probability. independent phases ipi a uniformly distributed in [0, 2tt] . 
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FIG. 8 Eigenvalue density of TV x TV Euclidean matrix A with elements Aij = exp(— |rj — rj|/f) — 8ij Ylk=i ex P( — l r * — r *l/f ) a ^ 
densities p£ 3 = 10~ 3 (left panel) and 10 -2 (right panel) for N = 10 3 points r< randomly distributed inside a cube. The analytic 
result ( 178 I (red dashed line) is compared to the numerical diagonalization (blue solid lines) and to the RG approach (green 
squares). The insets show p[ln(— A/2)] to emphasize the deviation from p(A) ~ 1/|A| that would correspond to p[ln(— A/2)] = 
const. 



account in Sec. II. F. 3 and corresponding to elementary 
excitations that involve clusters of many points. Unfortu- 
nately, the results of Sec. |II.F.3| cannot be used as such, 
since they were obtained for ERMs defined by Eq. ([9| 
with u = instead of u — 1. 



In addition to Eq. (76), p(A) of the exponential ERM 
can be studied with the help of the real-space RG ap- 
proach presented in Sec. |II.H| There, we argued that RG 
gives quantitative results in the low density limit p(, d <C 1 
only. This is confirmed in Fig. [8] where the symbols repre- 
sent the results of the RG flow following from Eqs. ( 167 1, 
(1681 and (169). These results are in agreement with 



those of Amir et al. (|2010|) who used a renormalization 

At low den- 



1671 



rule that is slightly different from Eq. 
sities p£ 3 < 1(T 3 , when both Eq. ( fl78[ ) and the RG pre- 
diction are good estimates of the eigenvalue density, the 
typical size m c (A) of the eigenvector associated with the 
eigenvalue A is, according to Eq . (|171 ), given by th e in- 
verse of the cumulative of (178) 



mir et al. 2010): 



m c (A) 



a 2 P e i [ln(-A/2)] 3 /3 



(179) 



This indicates that the eigenvectors are localized with a 
spatial extent m c (A) that diverges when A goes to zero. 
By means of a more refined analysis based on the percola- 
tion theory, [Amir et aT. (2012 1 show that a small fraction 
{p£?) v of delocalized eigenvectors survive close to A = 0, 
with v the critical exponent of the percolation transition. 



3. Cardinal sine Euclidean random matrix 

A common feature of the Gaussian and exponential 
ERMs considered in the two previous subsections is the 
fast decay of their generating function / (r i7 iv, ) with |r< — 
As a consequence, at a fixed density p of points rj and 
in the limit of large AT, the density of eigenvalues p(A) 
becomes independent of N and is controlled by a single 
parameter p. In this section we consider an example of 
ERM generated by a function / that decays slowly and 



oscillates in space, so that p(A) always depends on two 
parameters (e.g., p and N), even in the limit of N — > oo. 
Consider a matrix S with elements 



Sij = f(Ti ~rj) = 



sin(fco|r a 



(180) 



with a wavenumber fco = 2t:/Xq and a wavelength An (this 
terminology anticipates the application of the matrix S to 



and ' 


V.D.3 


. The matrix S was studied bj 


■ Akkermans 


et al. 


(2008) 


; Scully 


(2009|;|Svidzinsky et al. 


2008, 


2010); 



Svidzinsky and Scully (2009) in the context of light scat- 
tering and emission by large ensembles of atoms. 

An approximate solution for the eigenvalue density of 
the ERM ( 180 ) for points randomly distributed in a 



box of side L can be found by calculating the matrix 
T in the limit of large kgL using Eq. (13) and then ap- 



plying the free probability theory as d escribed in Sec. 
II.G.3| ( |Skipetrov and Goetschy| |2011| ). T can be ap- 
proximated by a diagonal matrix: T ~ (N/M)Im with 
M = (k L) 2 /2.8. This readily leads to the Marchenko- 
Pastur law for the eigenvalues A of the matrix S: 



p(A) = 



1 

1 

7 



5(A) 



V(A + -A)+(A-A_)+ 
2-7T7A 



(181) 

where A± = (1 ± -^7) 2 . This distribution is parameter- 
ized by a single parameter 7 = (A 2 ) — 1 ~ 2.8N/(koL) 2 
equal to its variance. It is easy to show that the same re- 
sult holds for arbitrary shape of the volume V, provided 
that 7 is calculated accordingly. In Fig. [9] we compare 
Eq. (181) with numerical simulations in a sphere of ra- 
dius R, where 7 = 9N/8(koR) 2 . T he ag reement is good 
at low densities and 7 < 1, but Eq. ( 181 ) fails to describe 



numerical results at high densities and 7 > 1. A better 
solution for p(A) is therefore required. 

Let us calculate the eigenvalue density of the matrix 
S using the self-consistent equation ( 125 1 derived in Sec. 



II. F. 3 As explained in Sec. II.F.3.d a general way to 
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FIG. 9 Eigenvalue density of the N x N ERM 5* denned by Eq. (180 1, where the N points r; a re ra ndomly chosen inside a 

rchenko 
27r/fc ) 



sphere of ra dius R. Numerical results (blue solid line) are compared to the Marchenko-Pastur law ( 181 1 (red dashed lines), and 
to Eq. ( 185 I (SC, green solid lines) at 4 different densities pAo of points (Ao 



solve this equation is to express the latter in the eigen- 
basis of the operator T. In order to solve the eigenvalue 
equation (135), it is convenient to decompose its kernel 



in spherical harmonics (Gradshteyn and Ryzhik 19801: 
sin(fco|r — r' 



fc |r 



= 47r X! jl( k or)ji(k r') 

1=0 m=-l 

xY lm (6,<t>)Y lm (9',<t>')*, 



(182) 



where 9 and <j) are the polar and azimuthal angles of the 
vector r, respectively, ji are spherical Bessel functions of 
the first kind, and Yi m are spherical harmonics. Inserting 
this decomposition into Eq. (135), we readily find that 



K a (r) = 1Z lm (r) = Am{k o r)Y lm (0,<t>), 
Ha = Hi = 4tt/? / dr' ' ji{k Q r') 2 r' 2 



(183) 



= 2^ LMfcoi?) 2 - ji-i(k R)ji+i(koR)] , (184) 

where Ai are normalization coefficients and a = {l,m}. 
Eigenvalues ^/ are (2Z + l)-times degenerate (m £ [—1, 1]). 
Eq. ( 136 ) then becomes 



1 



g(z) N 



V ^ 1 



(2i + l)w 



(185) 



Once this equation is solved for g(z), p(A) can be ob- 
tained from Eq. (l2lj). In contrast to Eq. (|181|), p(A) 



following from Eq. ( 185 ) depends on two parameters 7 



and pAq. 

As can be seen from Fig. [9j p(A) following from Eq. 
p5 l is in good agreement with numerical results at all 



densities pAg and for all values of 7. It even describes 
the splitting of the support of p(A) is separate 'islands' 
for R < Ao, even though the agreement with numerics is 
slightly worse in this regime, probably due to the finite 
value of N in numerical calculations [the limit N — > 00 



is taken in the analytic result ( 185 1] 



4. Cardinal cosine Euclidean random matrix 

Another example of ERM generated by a slowly de- 
caying and oscillating function f(vi,Yj) is the matrix C 
with elements 



C l3 = /(r, - Tj ) = (1 - Sij) 



cos(fc |r t — Tjj) 
fcn|r, : -rd 



(186) 



This matrix is relevant, for example, for understanding 
the collective Lamb shift in the problem of light inter- 
action with clouds of cold atoms (see Sec. IV.D.l). In 
contrast to Eq. (180), the Fourier transform of Eq. (186) 



is not always positive, and hence, the spectrum of C is 
not bounded from below. 

In a complete analogy with Sec. |II.I.3[ we compute the 
matrix T [see Eq. Q] for the ERM C and then apply 



the free probability theory of Sec. II. G. 3 (Skipetrov and 
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FIG. 10 Ei 


jenvalue density of the N x N ERM C defined by Eq. (|186|), where the N points Vi are randomly ch 


asen inside a 


sphere of radius R. Numerical results (blue solid lines) arc compared to Eq. (|187|) (red dashed lines), and to E< 


l. (|192p (SC, 


green solid '. 


ines) at 4 different densities pAg of points (Ao = 2n/ko). 











Goetschy 20111. This time, however, the calculation can 
be performed analytically only for the 7?.-transform of the 
probability density: 



is easy to show that the eigenvectors of f are necessarily 
of the form 



^/ s 2 , 47I " 3 7 

TZ(z) = arccoth — -= 

7T pAg 



1 



arctan 

2 



27T 3 7 



-1 - 



2tt 2 ' 



(187) 



arctan 



1 



PK 
2tt 2 ' 



1 



2tt 2 ' 



TZ a {r) = Kimpir) = Ai p ji(Ki p r)Yi m (6 ,<f>), 
where the coefficients Ki p obey 

ji(nipR) n ; _i(fc i?) 



k 



ji-i(nipR) ni(k a R) 



(189) 



(190) 



The resolvent g{z) has to be found numerically by invert- 
ing this equation. 

As can be seen from Fig. 



Integer p labels the different solutions of this equation for 
a given I. either real or imaginary numbers, and 

the corresponding eigenvalues 



10 Eq. (187) is not sufficient 



at high densities and a better solution is needed. We will 
search for it using the self-consistent equation ( 125| of 



Sec. II. F. 3 The eigenvalue equation (1351 is solved by 



Mo 



1 



2tt 2 



(Klp/ko) 



(191) 



using the expansion ( jGradshteyn and Ryzhik| |1980[ ) 
cos(fc |r - r'|) 



OO I 



fcolr 



= ~ 47r E E ji[komin(r,r')\ (188) 

;=0 m=— I 

x m [k m^{ry)]Y lm {6,4>)Y lm {e\4>')*, 



are (21 + l)-times degenerate (m g [—1 , 1]). In terms of 
the solutions m P of Eqs. ( 190 ) and ( 191 1, Eq. ( 136 ) reads 
finally 



1 



(192) 



where 9 and 4> are the polar and azimuthal angles of the 
vector r, respectively, Yi m are spherical harmonics, and ji 
and ni are spherical Bessel functions of the first and sec- 
ond kind, respectively. Inserting Eq. ( 188 ) into Eq. ( 135 ), 



We solve Eq. (190) numerically for K; p , then compute 
flip using Eq. (191), solve Eq. (192) numerically to find 



g(z) and, finally, obtain p(A) from the imaginary part of 
g(z) using Eq. (21 ). The results are shown in Fig. 



and using standard properties of spherical harmonics and 
spherical Bessel functions ( Morse and Feshbach 1953 ) , it 



see that in contrast to Eq. (187), Eq. (192) describes the 



10 We 



eigenvalue distribution quite nicely even at high densities, 
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but both equation fail for negative eigenvalues close to 
zero, where they predict a gap in the eigenvalue density, 
whereas numerical results do not show evidence for such a 
gap. Detailed analysis reveals that states corresponding 
to the eigenvalues in the gap are localized in space, sug- 
gesting a possible link between the failure of Eqs. (187) 
and (192| and Anderson localization (IGoetschy 2011). 



However, no definite conclusion can be made about this 
issue at the time of writing. 



III. NON-HERMITIAN EUCLIDEAN RANDOM MATRIX 
THEORY 



arbitrary non-Hermitian ERM in the limit of large ma- 
trix size (N — > oo) is developed in Sec. III.B Alternative 
approaches are also briefly discussed (Sec. III.C). We il- 
lustrate the theory by applying it to two specific non- 
Hermitian ERMs in Sec. IHI.DI 



A. Foundations of the non-Hermitian random matrix theory 

This section is devoted to the introduction of basic defi- 
nitions and relations useful in the study of non-Hermitian 
matrices. 



In contrast to Hermitian matrices, the eigenvalues of 
non-Hermitian matrices are not constrained to lie on the 
real axis and may invade the complex plane. Conse- 
quently, various methods developed for Hermitian ma- 
trices and heavily exploiting the analytic function theory 
are no longer applicable and require non-trivial modifica- 



tions 


Feinberg 


2006 Feinberg and Zee 


et al. 


1997a|b 


Jarosz and Nowak 2006) 



Most of the 



iterature on random non-Hermitian ma- 
trices focus on Gaussian randomness. A paradigmatic 
example is the ensemble of N x N matrices A with the 
probability distribution 



P(A) = C N e 



-NTrAA 1 



(193) 



Ginibre| (1965) showed that in the limit N — > oo, the 
eigenvalues of A are uniformly distributed within a disk 
of radius unity on the complex plane. Twenty years later, 
Girko ( 1985 ) generalized this result to matrix elements 
Aij that are i.i.d. with zero mean and variance 1/N. This 
is commonly referred to as Girko's law. Here we would 
like to tackle the problem of computing the density of 
eigenvalues of matrices that break away from this law: 
the non-Hermitian ERMs. Non-Hermitian ERMs appear 
in such important physical problems as Anderson local- 



ization of light ( Pinhciro et al. 2004 Rusek et al. 2000a ) 



and matter waves (Antezza et al. 2010 Massignan and 



Castin[ |2006"1), random lasing (|Goetschy and Skipetrov 



2011a||Pinheiro and Sampaio||2006[), propagation of light 



in nonlinear disordered media (Gremaud and Wellens 



2010), and collective spontaneous emission of atomic sys- 



terns (Akkermans et al, 2008; Ernst, 1968 


Skipetrov and 


Goetschy 201 1| Svidzinsky et al. |2010 


). However, no 



analytic theory was available to deal with these matri- 
ces until very recently, and our knowledge about their 
statistical properties was based exclusively on large-scale 



numerical simulations (jAntezza et al. 2010 


Gremaud 


and Wellens 


2010| Massignan and Castin[ 


2006 


Pinheiro 


et al.\ 2004| 


Pinheiro and Sampaio[ 2006 


Rusek et al.\ 


2000a Skipetrov and Goetschy| |2011|). 



Our review of non-Hermitian ERMs is organized as 
follows. In Sec. IIII.AI we introduce new mathematical 
objects that allow to generalize the methods developed 
for Hermitian matrices to the non-Hermitian case. A di- 
agrammatic theory for the density of eigenvalues of an 



1. Eigenvalue density and Hermitization 

Eigenvalues A„ of a N x N non-Hermitian matrix A 
are, in general, complex. Their density on the complex 
plane is 



p(A) = ^(f> (2) (A-AJ 



(194) 



where we use a shorthand notation <5^ 2 '(A — A„) = 
<5(ReA — ReA ?l )<5(ImA — ImA„). The relation between 
p{A) and the resolvent 



l(z) = 4 /ti 1 



-it — 



N \ z-A N \ ^— ' z — A r 

\n=l r ' 



(195) 



can be found using d z *(l/z) = Tr5(x)S(y), with the stan- 
dard notation d z * = \{d x + id y ) for z = x + iy. We 
obtain: 



p(A) - -d,.g{z) 

7T 



1 
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[d x Reg(z) - d y lmg(z) 



(196) 
(197) 



z=A 



Note that d y Keg(z) = —d x lmg(z) because p(A) is real. 
The right-hand side of Eq. (197) vanishes if g(z) obeys 
the Cauchy-Riemann conditions, i.e. if it is an analytic 
function of the complex variable z. In general, the eigen- 
values A„ occupy, on average, a two-dimensional domain 
T> on the complex plane where g{z) is nonanalytic, and 
p(A) describes the location and the amount of this non- 
analyticity. 

We now recall that the resolvent g(z) can be inter- 
preted as the electric field g(z) created at a point z on 
the complex plane by charges q — 1 situated at po- 
sitions An, see Eq. (44). Equation (197) can thus be 
seen as the Gauss law p(z) = ^7 x ,y • g(z)/2n. Hence, 
we readily obtain a new relation between p(A) and the 
logarithmic pairwise re puls ion V lnt (z) defined by g(z) = 



-V x , y V int {z) [see Eq. (143} 



p(A) = A xv V int (z) 

Fy ' 2irN ' y V ' 



(198) 
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where A. xy — Ad z d z - is the Laplacian in the coordinates 



x and y. Clearly, Eq. ( 198 1 may be particularly use- 



ful within the framework of the Dyson gas model where 
V ult may be related to the one-body potential determined 
by the probability distribution P(A) (see Sec. II. C and 
III.C.ll). 



In sertin g the explicit expression (43) of V lnt (z) into 
Eq. (198), we express p(A) in an alternative form: 



d z d z , (\ndet{z - A)(z* - 



1 

]™^<9 Z * (Indet [H A (z) ~ ithn] 



(199) 
(200) 



z=A 



(201) 



where I2N is the 2N x 2N identity matrix and Ha is the 
2N x 2N chiral Hcrmitian matrix 



(202) 



Note that Eq. ( 199 1 can also be derived from Eq. 
using d z d z *hizz* = irS(x)S(y). The representation 
is generally used in field-theoretical approaches 



194) 



201 



(Se^ 

III.C.2I). In addition, since the matrix (|202|) is Hcrmi- 



tian, one can compute its resolvent with well-established 
Hermitian techniques, from which it is still possible to 
recover the eigenvalue density of A ( |Feinberg[ 2006 Fein- 
berg and Zee| |1997b[ ). This is the so-called 'Hermitiza- 
tion method'. In the following, we will use an alternative 
method which has various advantages: it is technically 
slightly simpler, it reveals a relation between g(z) and 
the correlator of right and left eigenvectors of A, and, fi- 
nally, it allows for a generalization of the free probability 
calculus. 



2. Quaternions and the eigenvector correlator 

If A is Hermitian, its eigenvalues A n lie, on average, 
on some intervals (cuts) of the real axis. Therefore, it 
is possible to reconstruct g(z) by analytic continuation 
of its series expansion ( |27[ ) performed in the vicinity of 
\z\ — > 00. The eigenvalue distribution p(A) follows from 



the discontinuities of g(z) on the real axis [see Eqs. (21) 
and ( 196 )]. For a non-Hermitian matrix A, however, g{z) 
loses its analyticity inside a two-dimensional domain T> 
where A„ are concentrated, meaning that g{z) for z G T> 
cannot be simply assessed by analytic continuation of its 
series expansion. A way to circumvent this problem is 
based on the algebra of quaternions: while p(A) for an 
Hermitian A is obtained by approaching the real axis 
from orthogonal directions (in the complex plane), p(A) 
for a non-Hermitian A can be found by approaching two 
sides of V from directions 'orthogonal' to the complex 
plane in the quaternion space (Jarosz and Nowak 2006). 



Doubling the size of the matrix under study, we now work 
with a new 2N x 2A~ matrix 



a d_ A 
A - \ At 

and a quaternion resolvent matrix 



N\ Q®I N -A D 



(203) 



(204) 



The 2x2 matrix Q is an arbitrary quaternion in the 
matrix representation: 



x h + ix • a, 



(205) 



where x = (xi, X2, X3), <r is the triplet of Pauli matrices, 
a = xq + i^3, and b = x\ +1x2- Tr^ in Eq. (204) denotes 
the block trace of an arbitrary 27V x 2A^ matrix X. It 
is defined by separating X in four N x N blocks Xn, 
A12, X21, X22 and taking the trace of each of the latter 
separately: 



Tr N X = Tr 



An A12 
A21 A22 

TrAn TrAi 2 
TrA2i TrA22 



(206) 



Algebraic properties of the quaternions are useful to gen- 
eralize the free probability theory to non-Hermitan ma- 



trices (see Sec. III.C.3). However, if we wish to compute 
g(z) by a diagrammatic approach, it is sufficient to con- 
sider the quaternion Q — Z e , where 



z e = ^ j . (207) 

The generalized resolvent matrix G{Z t ) is then safely 



equal to its series expansion in 1/Z e (Janik et al. 1997b 



2001 Jarosz and Nowak 2006). By evaluating the block 



trace in Eq. ( 204 ) explicitly, one readily finds that 



with 

Gil 

^12 

so that 



1 

N 



G(Z e ) 



Tr 



11 ^12 

12 Wl 



(z- A)(z* - 4t) + , 



N 



Tr 



1 



(z - A)(z* - At) + ■ 



lim G(Z e ) = 

e->-0 



g(z) c{z) 
c{z) g{z)* 



(208) 

(209) 
(210) 

(211) 



Interestingly, the off-diagonal elements c(z) = 
lim e ^o G\ 2 yield the correlator of right \R n ) and 



2G 



left \L n ) eigenvectors of A (Chalker and Mehlig 1998 



Janik et al. 1999 1 : 



/ N 

C(Z) = -j- (j2( L n\Ln)(Rn\Rn)S {2) (z - A* 
\ri=l 

= Nc(z) 2 . 



(212) 



This shows that c(z) must vanish on the boundary ST> 
of the support of the eigenvalue density T>. In order to 
obtain p(A), one should compute G(Z e ) at finite e e K 
(by a diagrammatic or any other approach) , then take the 
limit e — > to extract g(z) from the diagonal elements of 
(pill, and finally apply Eq. (fl96 1. 



3. Biorthogonal basis of left and right eigenvectors 

For the sake of completeness, we recall here basic prop- 
erties of right \R n ) and left \L n ) eigenvectors of a non- 
Hermitian matrix A (or operator A). By definition, 



A\R n ) = A n \R„ 
(L n \A = A n (L ri 



Ai\L n )=A* n \L n ), 



(213) 
(214) 



meaning that \L n ) are the right eigenvectors of A^ . Ob- 
viously, A and A^ have complex conjugated eigenval- 
ues for det{A - A n I n ) = = det(^ f - A*J n ). Be- 
sides, |L„) and \R m ) are necessarily orthogonal because 
(L n \A\R m ) = A n (L n \R m ) = A m (L n \R m ). Assuming 
that the eigenvalues A„ are not degenerate, we normalize 
\R n ) and \L„) such that 



N 



(L n \R ri 



Et i* pi 



(215) 



Note that (R n \R m ) ^ S nm . Finally, the following prop- 
erties hold: 

I N = ^2\R n )(L n \ =£l L n><i2»|. (216) 

n n 

TrX = J2(Ln\X\R n ), (217) 

n 

where X is an arbitrary matrix. 



B. Diagrammatic approach for non-Hermitian Euclidean 
random matrices 



(a) 

H = JL 

rf- JL 



(b) 



MfJL 



\ItJL 




FIG. 11 (a) Diagrammatic representations of the matrices H , 
H\ A = HTH\ and A D . Full and dashed lines propagate 
in the bases {ri} and {ip a }, respectively (see Sec. II. A. 3 1; 
T = pA. (b) Diagrammatic notation for pairwise contractions 
Q and loop diagrams for any scalar x in the basis {i"i}, and 
for any operator X in an arbitrary basis {ip a }- 



statistics for the elements of H simplifies diagrammatic 
calculations but is not essential, contrary to the assump- 
tion of independence of different elements that may limit 
the validity of presented results. 



1. Derivation of self-consistent equations 

We start by expanding the 2x2 resolvent matrix G(Z e ) 
defined by Eqs. (204 1 and (207) in series in 1/Z e = 
(l/Z e )®I N : 



N 



Tr 



JV 



1 



— A D — 



(218) 



Inasmuch as H ia are i.i.d. Gaussian random variables, 
the result of averaging (...) over the ensemble of matrices 
H can be expressed through pairwise contractions ^ 
only. Diagrammatic notations, already introduced in Sec. 
|II.F.3| to evaluate efficiently the weight of different terms 
arising in the calculation, are reproduced in Fig. |ll| ^a) 
for clarity. The 'propagator' 1/Z e will be depicted by 



1 

z: 



>i 2 



M 2 

z* 



112 
T 2 2 



(219) 



Since each contraction ^ brings a factor 1/N, and 
each loop corresponding to taking the trace of a ma- 
trix brings a factor N [see Fig. [TT|b)] , only the planar 
rainbow-like diagrams, that contain as many loops as 
contractions, survive in the limit N — ¥ oo. Such dia- 



grams appear, for example, in Fig. 12 where we show 
the beginning of the expansion of the two independent 
elements of G(Z £ ) defined by Eq. ( [208 1. 

In the standard way, rather than summing the dia- 
grams for the resolvent, we introduce the 2x2 self-energy 



Our goal is to derive equations for the resolvent g(z) 
and the correlator c(z) of an arbitrary N x N non- 



Hermitian ERM A with elements Ai 



/( r i> r i) 



rj|A| r 7 ) in the limit of N — > oo. For this purpose, we 



follow Goetschy and Skipetrov (2011b) and mak e use of 



the representation A = HTH^ introduced in Sec. II. A. 3 
with the assumption that H has i.i.d. complex Gaussian 
entries satisfying Eq. The assumption of Gaussian 



Gf 1= _ + 
11 ii 



G 12 ~ 





FIG. 12 Diagrammatic expansion of the two independent 
elements of the matrix G(Z e ). 
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12 2 2 



1 



FIG. 13 Diagrammatic expansion of the two independent elements of the self-energy E(z7 e ). Braces with arrows denote parts 
of diagrams that are beginning of diagrammatic expansions of the quantities which the arrows point to. 




± n = f + ± n .Jfg$\...f + ± n .J£c$\...i' 



FIG. 14 The elements Efj and Ef 2 of the matrix E(Z e ) 
can be written as traces of operators E^ and E 12 : E n = 
TrEfi/JV and Ef 2 = TrE E 12 /iV. Operators En = lim £ _ >0+ E xl 
and Ei 2 = lim (! ^. + E 12 obey coupled equations, where g = 
lim e ^ 0+ Gli and c = lim e ^ 0+ Gf 2 [see Eq. \211\]. 



matrix 



E(Z 6 ) = Z e - G(Z e 



12 

■^12 ^11 



(220) 



It is equal to the sum of all one-particle irreducible dia- 
grams contained in 



Z e G(Z e )Z £ 



1 

A 



Tr 



N 



A 



D 



A D —A D 
Z f 



(221) 

The first dominant terms that appear in the expansion 
of the two matrix elements E^ and £f 2 are represented 
in Fig. [13] In the two series of Fig. [13] we recognize, 
under a pairwise contraction, the matrix elements Gf-L 



as well as the two operators 
14 Equations obeyed by the 



and G\ 2 depicted in Fig. [12 
T,\-i and £f 2 denned in Fig. 
operators En — lim e _j. + £f i and £12 = lim^o £f 2 are 
obtained after summati on of all planar rainbow diagrams 
in the expansion of Fig. 13 and taking the limit e — » + . 18 
The diagrammatic representation of these equations is 
shown in Fig. |14| Applying 'Feynman' rules defined in 



Fig. [TTj^b) , we obtain: 



£11 = (l + .g£n + c£i 2 )T, 
£ 12 = (c£ 11 + 5 *£ 12 )ft ) 



(222) 
(223) 



where T = pA. After some algebra, 19 £n = TrEn/A 
and £12 = Tr£i2/iV can be expressed as: 



£11 



£12 — 



i Tr (i-g*rt)f 

A (l- 5 *ft)(l- 5 T) -c 2 Ttf' 

Tr ^ 
N 1 (l- ff *T , t)(l-gf)-c 2 ftf' 



(224) 
(225) 



Furthermore, as follows from Eq. (211 1 and the defini 



tion ( 220 1 of the self-energy matrix, g and c are simply 



related to En and £12 by 



g(z) c(z) 
c{z) g(z)* 



z -En -£12 
-£12 7* - Ej! 



(226) 



Elimination of the self-energy from Eqs. (224 1, (225) and 
( 226J) leads to two self-consistent equations for the resol- 
vent g(z) and the eigenvector correlator c(z): 



+ TrTr 



(l-g*f*)f 



C 2 A (l-g*ft)(l- 5 f)-C 2 ftf ' 



= ^Tr- 



A (l- 3 *Tt)(l- 5 T)-c 2 TtT 



(227) 
(228) 



As usual in such a procedure, summation must be performed 
before taking the limit e — > 0. Hence, the off-diagonal element 
of the propagator 1/Z e gives rise to non- vanishing terms after 
summation, although it is zero in the limit e — > 0. 



19 Although [T,Tt] 7^ 0, this calculation is easily performed by 
making cyclic permutations of operators under the trace opera- 
tor. 
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At this final stage, it is convenient to define the follow- 
ing operators 



So = S(g) = 



T 



i-aT' 



(229) 



Si = S(g + c 2 S^ 



(1 - g*T*)f 



(l-g*Tt)(f- 5 T)-c 2 TtT 

(230) 



in terms of which Eqs. (2271 and (2281 become 



* 1 

1 1 

I-/ - - i- N 



1 _ ,.2 _ i\r TrS ' 1 ' S 'o- 



(231) 
(232) 



Because c(z) must vanish on the boundary ST> of the 
support of the eigenvalue density T>, equations for z 6 ¥D 
follow: 



2. Analysis of self-consistent equations 



(233) 
(234) 



Equations (227 1, (2281, (233) and (234) are the main 



results of this section. An equation for the borderline of 
the support of the eigenvalue density of a non-Hermitian 
ERM on the complex plane z = A follows from Eqs. ( 233 ) 



and (234) upon elimination of g. The density of eigenval- 



ues A inside its support T> can be found by solving Eqs. 
(227) and (228) with respect to g(z) and then applying 
Eq. (1961). 



At low density and in the framework of representa- 
tion A = HTH\ the eigenvalue density of any Hcrmi- 
tian ERM obeys the Wigner semicircle law (131 1. 20 An 



analogous result exists for non-Hermitian ERMs as well. 
Indeed, using the approximation 



Si — So — T 



(235) 



valid at low densities, Eqs. (233) and (234) for the bor- 



derline of the eigenvalue domain reduce to 

2 



1 

AT 



-TrT 



= — TrTT t , 
N 



(236) 



become 

9{z) 



J_ TVJ't 



i TrTTt 



<z) = — 



1 



N 



TrTTt 



k-^Trf| 2 
i Trfft 



(237) 
(238) 



The term TrT/ A that appears in Eqs. |236j, (|237|, and 
( 238 1 , leads to a shift of the eigenvalue distribution equal 



to 



TrT 



TrA 



(Tr N A) 
N 



(A). 



(239) 



We assume from here on that An = (i = 1, . . . , N), so 
that, in particular, (A) = 0. With this assumption, the 
term Trfft rea d s: 



Tr(f f t) = p 2 Tr(ii f ) = 



d d r dV|/(r,r')r 



(Tr N (AA*)) 

INN \ 

( E J2 KKn(L n \L m ){R m \R n ) 



\n— 1 m— 1 
/ N 



(J2\ A n\ 2 (L n \L n )(R n \R n ) 

\n=l I 

2/ElA„| 2 \=2A(|A| 2 ). 



(240) 

(241) 
(242) 
(243) 



\n=l 



In Eqs. (242) and (243) we assumed that, at low den- 



sities, (rj|L ra ) and (ri\R n ) behave as Gaussian random 
variables. Introducing the shorthand notation 



1 = 



Tr(fft) 



2A 



<|A| 5 



we rewrite Eqs. ([236]) , ( |237[ ), and ([238]) as 

\z\ 2 = 2 1 (ze5V), 



9(z) = J- (« e 23), 



c{z) 



2" 
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1 



(zev). 



(244) 

(245) 
(246) 

(247) 



This shows that, in the limit N — > oo and p — > at 
fixed 7, the eigenvalues of an arbitrary traceless non- 
Hermitian ERM are uniformly distributed within a disk 
of radius v / 2~7- Within the disk, p(A) = 1/2ttj. This is 
the famous Girko's law (Girko 1985), first discovered by 



and Eqs. ( |231[ ) and ( |232| ) for g(z) and c(z) with z e D |Ginibre| ( |1965[ ) for the complex Gaussian ensemble. We 



20 We exclude rare ERMs for which the operator T has a small 
number of non-zero eigenvalues. 



recover this law because in the limit N — > oo and p — > 0, 
elements of A essentially behave as i.i.d. variables and 
hence En = and E12 = c(z) /2j. 

As it was the case for Hermitian ERMs, the solution 
of Eqs. p27) , p28| , p33| and p34| for a given matrix 
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A can be greatly facilitated by a suitable choice of the 
basis in which traces appearing in these equations are 
expressed. In addition to {r} and {k Q }, a biorthogonal 
basis of right \TZ a ) and left \C a ) eigenvectors of T can 
be quite convenient. We recall that the right eigenvector 
\K a ) obeys 

(r\f\K a ) = p [ dV/(r,r')K a (r') = fi a TZ a (r), (248) 
Jv 

where /i Q is the eigenvalue corresponding to th e eigenvec- 
tor \lZ a )- The traces appearing in Eqs. (233) and (234) 
can be expressed as 



fJ>a 



g^c 



a,0 



(1 ~ 9Ha)Q-- 9Pp)* 



(249) 



(250) 



respectively. Technically, the main difference with the 
study of Hermitian ERMs is that we now hav e to k now 
the eigenvectors of T explicitly [compare Eqs. (250) and 
p6t]. 



C. Other approaches 

1. Dyson gas picture 

In this section, we extend the Dyson gas picture intro- 
duced for Hermitian matrices in Sec. |II.C[ to the non- 
Hermitian case. Let us first consider the ensemble of 
non-Hcrmitian random matrices with Gaussian probabil- 
ity dmsity defined by Eq. ( 193 ) and originally introduced 



by Ginibre (1965). The probability density P(A) of this 
ensemble is invariant under all unitary transformations, 
but not under the similarity transformation A —> SAS~ 1 
used to diagonalize A = S~ 1 DS (D denotes a diagonal 



matrix with elements A„ 



1...JV). Hence, P(A) de- 



pends explicitly on S and not only on the eigenvalues of 
A. This feature is the main difference with the Gaussian 
ensemble or the Wigner-Dyson ensemble ( |40| . Her- 
mitian matrices drawn from these two latter ensembles 
can be diagonalized by unitary matrices, so that P(A) 
depends on {A n } only. In order to obtain the joint prob- 
ability density P({A„}) from Eq. (193), we must change 



variables from Ay to parameters related to eigenvalues 
and eigenvectors of A. Since TrAA^ depends on eigen- 
vectors, the new variables have to be chosen carefully to 
facilitate further manipulations. The result is the follow- 
ing (|Ginibre| |l965{ |Mehta| [2004]) : 



P({A„}) = CV-^ 9({A "» 



(251) 



H3({\ n }) = Nj2v 9 (A n )- ^ln|A n -A m |, (252) 



V 9 (z) 



n=l 

2 ' 



(253) 



and /3 — 2. We recognize in Eq. (251) the Boltzmann- 
Gibbs distribution of a Coulomb gas in thermal equilib 

1/13. Equations (|25T]»-(|253|> 



rium at a temperature T = 1//3. Equations 
have exactly the same form as Eqs. (41 ) and (|42[). As for 
Hermitian matrices, the logarithmic pairwise repulsion 
comes from the Vandermonde-type Jacobian |V({A„})|°. 

In the limit N — > oo, we can perform coarse-graining 
of the energy functional H 9 [see Eq. ( 54 )] , and minimize 
it to obtain the equality 



d z V int (z) = Nd z V°(z), 



(254) 



where V (z) is the logarithmic pairwise repulsion (43) 



Eq. (254) means that the force Nd z V 9 (z) experienced by 
each particle of the gas is compensated by the Coulomb 
repulsion by all other particles. The eigenvalue distribu- 
tion ( 198 ) reads now: 



p(A) = ^z, v V(z) 



(255) 



z=A 



Note the difference with Eq. ( 53 1 for Hermitian matrices 



Eq. (255) is local; the shape of the distribution at A de- 



pends on the profile of V 9 in the vicinity of A only, while 



in Eq. (53) the shape of p(A) strongly depends on the 



boundaries of the distribution, meaning that the influ- 
ence of V 9 on p(A) is nonlocal. Note also that Eq. ( 255 ) 
contains no information about the borderline of the sup- 
port of eigenvalues. If V 9 is simple enough, the border- 
line can be obtained from the normalization constraint 
/ dAp(A) = l. 21 For V 9 (z) = \z\ 2 /2, we find that the 
eigenvalues are uniformly distributed inside a disk of ra- 
dius 1. This is the celebrated Ginibre 's result (Ginibre 



1965) 



Obviously, Eqs. (251) and (252) also apply to any nor- 
mal matrix A ([A, A 1 ] = 0) with probability density 



P(A) = C N e 



-NTrV a (AA 



(256) 



where V 9 is arbitrary, for A can be diagonalized by a 
unit ary m atrix. The one-body potential appearing in 
Eq. ( |252| is then given by V 9 ( z) = V 9 (\z 2 ). A counter- 
intuitive result is that solutions ( 251 ) and ( 252]) may com- 
pletely break down for most of random non-Hermitian 
matrices — i.e. for non-Hermitian matrices that arc not 



normal or partially normal (Feinberg 2006 Feinberg 



et at 



Fcin- 



2001 1 — distrib uted according to Eq. (256). 
berg and Zee ( 1997a[ ) proved the 'single-ring theorem'. It 
stipulates that the shape of the eigenvalue distribution is 
either a disk or an annulus, whatever polynomial the po- 
tential V 9 is. This is clearly in contradiction with what 
we could expect from Eqs. (12511) and 12521) , that tell us 



that the number of domains occupied by the eigenvalues 
on the complex plain should grow with the number of 



21 For co mpli cated cases, the borderline may be found by inspection 
of Eq. |55l that is still valid for non-Hermitian matrices. 
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minima of V g (z) = V 5 (|z| 2 ). The polynomial V 9 ~ AA^ 



that corresponds to the complex Gaussian ensemble ( 193 ) 



is actually the only polynomial for which Eqs. (251 ) and 
(252 1 are valid whatever the matrix A obeying (2561 is. 



Remarkably, Feinberg and Zee (1997a) also showed that given 



the eigenvalue distribution of A can nevertheless be found 
from the resolvent of the Hermitian matrix AA^ . This re- 
solvent has already been known in the literature for an 



arbitrary polynomial V 9 (Feinberg and Zee 1997a) 



Although the Dyson gas picture was not rigorously jus- 
tified for ERMs, it may be helpful to qualitatively un- 
derstand the eigenvalue distributions obtained by other 
methods. For example, in the study of the Green's matrix 
G in Sec. III.D.2 we observe that the support T> of p(A) 
deforms when the density is increased, going through a 
transition from a disk-like to an annulus-like shape, and 
eventually splitting into multiple disconnected domains 
at high density (see Fig. 16 1 . It is difficult to refrain from 
interpreting such transitions as phase transitions for the 
Dyson gas due to modifications of a hypothetic one-body 
potential V 9 . 



2. Field representation 

Let us now briefly explain how to compute the eigen- 
value distribution of a non-Hcrmitian matrix A in the 



field-theoretical approach. We start with Eq. (201) 
rewritten as 



p(A) 



— \im d z d z , (lnZ e (z)) 

TT1\ e->0 



z=A 



where we introduced the partition function 
Z e = det 



el n i(z - A) 
i{z* - At) el n 



(257) 



(258) 



In order to evaluate (\nZ e (z)), we follow the same proce- 
dure as in Sec. |II.E[ namely, we apply the replica trick, 



(lnZ e (z)) = lim 



together with the representation 



Z\z) OC J d(f>! 
N 



t> N e 



(259) 



(260) 



z, e) = ^ <t>\ {th + ixa x - vyu v ) 4>% 



i=l 
N 



(261) 



where the N fields 4>i are pairs of complex variables, a x 
and cjy are Pauli matrices, z = x + iy, and A = A h +iA s , 
with A h and A s Hermitian matrices. This representation 
combined w i th the cavity method was used by |Rogers| 
and Castillo ( 2009 ) to analyze the eigenvalue distribution 



of sparse non-Hermitian matrices. A slightly different 
representation of p(A), also based on the replica trick, 
can be found in a nice review of RMT by Ste phanovj 
et al. (2001), where a derivation of Girko's law is also 



For non-Hermitian ERMs of the form Aij — /(r, — r^ ) 
it seems feasible to generalize the field method proposed 
by |Mezard et al.\ for Hermitian ERMs. Basically, 

it amounts to make the same approximations in the cal- 
culation of (Z e (z) n ) as those presented in details in Sec. 
|II.E| This leads to equations that have the same degree 
of validity as Eq. (97). For example, the equations for 



the borderline of the eigenvalue domain are 



1 



d d k 



/o(k) 



(2tt)<* 1 - pf (k)g(z) 
Pl/ (k)| 2 



d d k 



(27r)*|l-p/o(k) fl (*)|* 



(262) 
(263) 



where /o(k) is the Fourier transform of / (r). These equa- 
tion can also be obtained from Eqs. (233) and (234) 
by using the approximation (kLA|k') ~ (k|A|k)5kk 
/o(k)5kk'- Contrary to Eqs. |233|) and p34|, Eqs. J262|) 



and ( 263 ) are parameterized on 
N/V. 



3. Free probability 



y by the density p = 



The extension of free probability theory and, in partic- 
ular, the generalization of the concept of Blue function, 
to non-Hermitian matrices is natural in quaternion space 
(Jarosz and Nowak 2004). The quaternion Blue matrix 
j 3x{Q) of any random matrix X is the functional inverse 
of the quaternion resolvent matrix ( 204[ ): 



G X [B X (Q)] = B X [G X (Q)] = Q, 



(264) 



where Q is a quaternion defined by Eq. ( 205 1 . For con 



venience, we also introduce the quaternion i?-transform: 

1 



R X {Q)=B X {Q)- 



Q' 



(265) 



For Q = Z e given by Eq. (|207|, R x {Z t ) is simply related 
to the self-energy matrix (220) by 



R x (Z e ) = -E x [B x (Z e )], 



(266) 



and therefore R x (z) = lim € _>.o R x (Z e ) and T, x (z) 
lim £ ^o ^x(Z e ) are related through 22 



R x {z) = Y, x [B x {z)], 



(267) 



where B x {z) is the usual Blue function (29). We now 
mention two important properties of the matrices G x (Q) 



To obtain Eq. ( 267 1, we use B[diag(^, z*)] = diag[B(«), S(z*)l. 
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andi?x(Q)- Firs t, Gx(Q) and Rx(Q) obey the fo llowing 
scaling relations (Jarosz and Nowak 2004 2006): 



G a x(Q) — Gx 
Rax(Q) - 



1/a 
1/a* 



a 
a* 



Rx 



Q 



Q 

a 
a* 



1/a 
1/a* 



(268) 
(269) 



where a G C*. Second, both Gx(Q) and Rx{Q) can 
be expressed in terms of the resolvent gx (z) and the 1 Z- 
transform TZx(z) only (Jarosz and Nowak{ 2004 2006): 



Gx(Q) = {[qgx(q) ~ q*gx(q*)] h 

~ i9x(q)-9x(q*)}Q f }, 

Rx(Q) = {[<Z^x(<7) - q*n x (q*)] h 

q-q* 

- [nx(q)~n x (q*)]Q^} 1 



(270) 



(271) 



where q = xq + i|x| and q* are two complex conjugated 
eigenvalues of Q. 23 

For arbitrary Q, we can use algebraic properties of 
quaternions to show that the following addition law holds 
(I Jarosz and Nowakl 120041 120061): 



R Xl+ x 2 (Q) = RxAQ) + Rx 2 (Q), 



(272) 



where X\ and X2 are two non-Hermitian, asymptotically 
free random matrices. Therefore, applying successively 
Eqs. (p9l, poll, p7il, p65j) and p64| for Q = Z e , we 



can infer Gx 1 +x 2 {Z7) from gx 1 {z) and gx 2 ( z )- The steps 
of the algorithm are: 



gx, -> S Xl -> K x , -4- - 



(273) 



The resolvent <?Xi+x 2 (z) and the eige nvec tor correlator 
c Xi+x 2 ( z ) are finally found from Eq. (211). 

This algorithm is greatly simplified when we look for 
the eigenvalue distribution of a non-Hermitian matrix 
Xi + 1X2, where X\ and X2 are free Hermitian matri- 



ces with known TvL-transforms. Jarosz and Nowak (2004 



2006 ) showed that the problem reduces to solving a sim- 



ple system of three equations with three unknown vari- 
ables, complex u, v, and real t: 



Kxi(u) =x + 

Kx2{v) = y- 
\u\ = \v\, 



t-1 



(274) 



The relation l |27l| l between Rx(Q) and TZx(g) holds also be- 
tween Bx(Q) and Bx{q) because Q _1 = Q< /qq* . 



where z — x + iy. We express u and v via t from the 
first two equations, substitute the results into the third 
equation, and then solve for t. The resolvent and the 
correlator are then given by 



gx 1 +ix 2 (z) = Re u - iRe v, 



(275) 



c Xl+i x 2 (z) = (Reu) 2 + (Rev) 2 - \u\ 2 . (276) 

Equation for the borderline z £ 8T> of the eigenvalue do- 
main follows from cx 1 +ix 2 {z) — 0. 24 From this simplified 
algorithm it is straightforward to recover the Girko's law 
for Gaussian Hermitian matrices X\ and X2 \TZ\(z) — 
T^i{z) = 2, see Eq. (98)]. A less trivial example of appli- 
cation of this algorithm is given in Sec. |III.D.1| 

Finally, a generalization of the concept of 5-transform 
( 147 1 to non-Hermitian matrices was presented by Burda 



et al. (2011). These authors derived the multiplication 



law analogous to Eq. ( 149 1 for Hermitian matrices, that 
allows to compute the eigenvalue distribution of a prod- 
uct X1X2 of two random non-Hermitian matrices X\ and 
X2 from the known properties of the latter. 



D. Application of the general theory to specific Euclidean 
random matrix ensembles 

1. Cardinal cosine +ix (cardinal sine) matrix 

We start our study of non-Hermitian ERMs by the case 
of a N x N matrix 



X 



= (1 - Sij) 



cos(fc |ri 



r,-|) .sin(fc |r£-r^. 



feo r< 



(277) 



where {r^} and {r^} are two different and independent 
sets of points. We recognize in the real and imaginary 
parts of X the two Hermitian ERMs C and S stud- 
ied independently in Sees. |II.I.4 and II. 1. 3 respectively. 
The matrix X = C + i[S' — In] is similar to the three- 
dimensional free-space Green's matrix G to be studied in 



Sec. III.D.2 except that its real and imaginary parts are 
not correlated. Using the definition ( 140 ) of asymptotic 



freeness, it is easy to check that the matrices C and S' 
are asymptotically free, in agreement with the intuitive 
definition of freeness as statistical independence. 

Since X is of the form X\ + 1X2, where X\ — C and 
X2 = S' — Jjv are two asymptotically free Hermitian ma- 
trices, we can make use of Eqs. (274), (275) and (276) to 



24 If we are interested only in the borderline 8T>, i.e. the bound- 
ary between the holomorphic and nonholomorphic domains of 
gx i+ix 2 (z), it is also possible to use a conformal transformation 
that maps the cuts i 6 1 of 9x 1 +x 2 (t) onto ST>. The equation 
z = f(t) follows from gx,+ix, (z) = 9x 1 +x 2 (t) ( |Janik et al\ 
1997a| Jarosz and Nowakl 20061. 
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FIG. 15 Densit y plo t of the logarithm of the probability density of eigenvalues A n of a square N x N Euclidean matrix X 
defined by Eq. (277 1 at 4 different densities p of points rj, per wavelength Ao = 2-7r/fco cube. 2N = 2 x 10 4 points 
and rj (i = 1, ...,N) are randomly chosen inside a 3D cube of side L; j — 2.8N/(koL) 2 . The probability distributions are 
estimated from 10 realizations of {r^} and {r£}. Dashed lines show the domain of existence of eigenvalues following from the 
free probability theory. [Reproduced from (Skipetrov and Goetschy 2011k] 



calculate the resolvent g(z) and the eigenvector correla- 
tor c(z) of X. In the limit of 7 <C 1, the ^-transforms of 
X\ and X2 are those of Gaussian and Wishart matrices, 
resp ectively : 72-Xi (z) = 7% and l Zx 2 (z ) = 1 /(1 — 72) (see 
Sec.|II.G.2|. Solving Eqs. p74|, ([275]) and p76|, we find: 



g(z = x + iy) 
c(z = x + iy) 



x 1 

2^ ~ 2 



1 





1 




" 4 



7(1 + y) 2 + y 

y 



, (278) 
2 



7(1 + 2/) 2 + y 



7(l + y)(2 + y)' 



(279) 



The correlator ( 279 ) must vanish on the borderline 6T> 



of the eigenvalue domain. We therefore readily obtain an 



equation for the borderline on the complex plane: 

2 



V 



7 



4t) 



= 0. (280) 



,1 + y 2 + y J (l + y)(2 + y) 
The probability density inside the domain delimited by 



Eq. (280) is 



p{x,y) 



1 

2~TX 

1 

4tt 



[d x Reg{z) - d y lmg(z)] 
"11 1 

_7 + 7(i + y) 2 " W+W. 

A better model for the "^-transform of the matrix 



(281) 



X\ = C is given by Eq. (1871. If we use this equa- 
tion instead of 1Zx 1 (z) = 72 above, analytic calcula- 
tion becomes impossible but we can still compute g(z) 
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and c{z) numerically. The resulting borderline of the 
eigenvalue domain is shown in Fig. 15 (dashed lines) 



together with the eigenvalue distribution of the matrix 
X = C + i[S' — In] found by the numerical diagonaliza- 
tion of a set of 10 4 x 10 4 random matrices. At the lowest 
density considered pAg = 0.01, the borderline following 



from Eq. (187) is very close to Eq. (280). At higher den- 
sities, the former describes numerical results much better 
than Eq. ((2801). 



Equation (280) predicts a splitting of the eigenvalue 
domain in two parts at 7 = 8. The more accurate calcula- 
tion using Eq. ( |187[ ) leads to a similar prediction (see the 
lower right panel of Fig. 15). However, the eigenvalues of 



the matrix X do not show such a splitting and form an 
'inverted T' distribution on the complex plane instead. 
This is due to the fact that the Marchenko-Pastur law 



( 181 1 fails to describe the eigenvalue distribution of the 
matrix S' at 7 > 1 and hence the 7?.-transform 1/(1 — 72) 
that we assumed for S' is not a good approximation any- 
more. 



2. Random Green's matrix 



Let us now illustrate the power of Eqs. (227), (228), 
(|233| , and ( 234 ) on the example of the N x N random 



Green's matrix 



Gij = (!-%) 



cxp(ifc |rj - r 3 |) 



(282) 



This non-Hermitian ERM is of special importance in 
the context of wave propagation in disordered media be- 
cause its elements are proportional to the Green's func- 
tion of Helmholtz equation, with that may be thought 
of as positions of point-like scattering centers. It ap- 



peared in works of Ernst (|1968); 


Goetschy and Skipetrov 


(|2011a|b|);|Gremaud and Wellens| Q2010D; |Massignan and 


Castin| (|2006|) ; Pinheiro et aL|(|2004|); Pinheiro and Sam- 


paio (2006); 


Rusek et al. 


2000a) 


Skipetrov and Goetschy 


(2011); Svic 


zinsky et al 


(2010 


), but was studied only 



by extensive numerical simulations, except in the pa- 
per by Svidzinsky et al. (2010) where analytic results 



were obtained in the infinite densi ty limit, and our works 
( Goetschy and Skipetrov 2011a|b ) where an analytic the- 
ory applicable at any density was developed. 

Very generally, the eigenvalue density of G depends 
on two dimensionless parameters: the number of points 
per wavelength cubed pAg and the second moment of 
|A| calculated in the l imit of low density: (|A| 2 ) = 7 = 
Tr(ff*)/N [Eq. |244| ], Even though the latter result 
for (|A| 2 ) can be rigourously j ustifi ed only in the limit of 
low density pAg <C 1 [see Eq. (243)], we checked numeri- 
cally that it holds approxi mately up to densit ies a s high 
as pAg - 100. Equations {243)), p44| and (p0| show 



moments of the eigenvalues of ReG and ImG: 



7 



(|A G | 2 ) = ((ReA G 



((ImA G ) 2 



( a Lg) 



(ALg) 



(A 2 



mGI- 



(283) 
(284) 



Eq. (2841 holds for k R > 1 and pAg < 100, whereas 
Eq. (283) is also valid for any k R and arbitrary non- 



Hermitian traceless ERM [provided that Eq. ( 243 ) holds] . 
We will see from the following that the two parameters 
pAg and 7 control different properties of the eigenvalue 
density. 



a. Borderline of the eigenvalue domain: approximate solution 
at low density. We first focus on the borderline of the 
support of eigenvalues which is easier to visualize. We 
assume that the N points are chosen inside a sphere of 
radius R. For arbitrary fcgi?, the parameter 7 is then 
given by 7 = 9N/8(kpR) 2 . Traces appearing in Eqs. 
(233) and ( 234[ ) in the (r)-representation read 



TrSg = Tr 




T 



= Tr [T + gTSi 
1-gTJ V 

d 3 r d 3 r'T(r, r')5 (r',r), 
Tr^ = ^d 3 rdV|S (r,r')| 2 , 



(285) 
(286) 



where T(r, r') = p(r|A|r') = pexp(ifc |r — r'|)/fc |r — r'| 
and in Eq. (jggs} we used the fact that Trf = pTri = 0, 



as follows from Eq. (|282j). S (r,r') = (r|S |r') obeys 

S (r,r')=T(r,r')+5 / d 3 r"T(r, r")S (v", r'), (287) 
Jv 

as follows from the definition of Sq. Noting that 

(A r + fc 2 + ie) T(r, r') = -^6^(r - r'), (288) 



where e — > + , we apply the operator A r 



ie to Eq. 



( |287| ) and obtain 
A r ^ (r,r') + fc 2 



So(r,r') 



47rp 



^(r-r'), 



(289) 



where lly(r) = 1 for r € V and elsewhere. In the 
limit of low density pA 3 , — > 0, an approximate solution of 
this equation is obtained by neglecting 'reflections' of the 
'wave' S'o(r,r / ) on the boundaries of the volume V and 
thus setting Ily(r) = 1 everywhere. This yields 

exp [ift(<7)|r — r'|] 
fcolr-r'l 



that the second moment of |A| is related to the second 



S (r,r') 



P- 

k \l I 



2tt 2 ' 



(290) 
(291) 
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FIG. 16 Density plots ol the logarithm of eigenvalue density of the N x N random Green's matrix G obtained by numerical 
diagonalization of 10 realizations of the matrix for N = 10 4 . Points i\ are randomly chosen inside a sphere of radius R. The 
solid red lines represent the borderlines of the support of eigenvalue density following from Eq. ( 296 I in panels (a) and (b) and 
from Eqs. (3141 and (315 I in panels (c) and (d). The dashed lines show the diffusion approximation (300 1. [Reproduced from 



(Goetschy and Skipetrovl 2011b I 



We now plug the explicit expressions for T(r, r') and 
S (r,r') into Eqs. ([285]) and dgggb . This yields 



Tr5 = 2jNgh[-m(g)R - ik R], 
TiSoSl = 2jNh[2ImK(g)R], 



(292) 
(293) 



with 

h(x) 



6x 4 



[3 - 6x 2 + 8x 3 - 3(1 + 2x)e- 2x ] . (294) 



In th e low- dens ity limit, g can be eliminated from Eqs. 
( p33) and |234| by neglecting Tr S /N in Eq. |233| and 
substituting g = 1/z into Eq. ( 293[ ). This gives 



|A| 2 = 2-fh [2Im K (l/A)i?] . 



(295) 



If the argument of th e fun ction h in Eq. ( 295 1 is expanded 
in series in pAp, Eq. (295) becomes: 



|A|^ ~ 2-yh -87 



ImA 

3|Af 



(296) 



By comparing Eq. ( 296 ) with the exact solution (see Sec. 
III.D.2.C and Fig. 16), wc conclude that it is valid up to 
densities as high as ~ 10. 

For 7 <C 1, the density of eigenvalues is roughly uni- 
form within a circular domain of radius 1/27, see Fig. 



16 



a). The domain grows in size and shifts up upon in- 
creasing 7. At 7 > 1 it starts to 'feel' the 'wall' ImA = — 1 
and deforms [Fig. 16 'b)]. Before considering the shape of 
the eigenvalue domain at higher densities, we would like 
to show how the link with scattering theory can be used 
to derive an equation for its borderline. 



b. Borderline of the eigenvalue domain from the scattering 
theory. Because the elements Gij of the Green's matrix 
G describe propagation of a scalar wave between points 
Ti and Tj in space, it is quite natural that a link exists 
between some of the properties of the eigenvalue den- 
sity of G and the theory of wave scattering in an en- 
semble of point scatterers. In particular, if we define 
Iij as the intensity of a wave at Vj due to a source at 
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Ti, I(t) = ^2i£j Ijj can be written as (Goetschy and 
Skipetrov[ [201 lb| ) 



J(t) = IV 



(297) 



where t is the scattering matrix of an individual scatterer 
and Q = —koG/Air. This is to be compared with the 
expression for the correlator of right and left eigenvectors 
of an arbitrary matrix A, c(z) = lim e _ i .Q+ G\ 2l following 
from Eq. pTo|: 



c(z) 



ie / 1 

e ™+iV\ r (z - A){z - A)t + e 2 



For A = Q 1 and z = t we thus have 



c(t) = 



(298) 



(299) 



This should become different from zero when i enters the 
support of the eigenvalue density of Q~ x or, equivalently, 
when 1/t enters the support of the eigenvalue density of 
Q . The only way to obtain c(t) ^ for e — >• + is to make 
{I(t)) diverge. In the framework of the linear model of 
scattering, this can be achieved by realizing a random 
laser (Goetschy and Skipetrov 2011a). We thus come 



to the conclusion that finding the borderline of the sup- 
port of the eigenvalue density p(A) of the N x N Green's 
matrix ( 282 ) is mathematically equivalent to calculating 
the random lasing threshold in an ensemble of N identical 
point-like scatterers with t = — Air/koA. Application of 
this link to the study of random lasers will be discussed 
in Sec. |IV.D.2| In the diffusion approximation, for ex- 
a mple, the threshold of such a random laser was found 
by Froufe-Perez et al. ( |2009 ) and leads to the following 
equation for the borderline: 



|A| 2 



87 



Vl + ImA 1 



|Ap 



|A| 2 +4 7 



(300) 



We show this equation in Figs. [16[a) and (b) by dashed 
lines. It gives satisfactory results only in the weak scat- 
tering regime pAg < 10. However, it should be noted 
that in this regime, it describes the lower part of the 
borderline, corresponding to small ImA, better than Eq. 
p96|. 



c. Borderline of the eigenvalue domain: exact solution of self- 
consistent equations. The approximate equation ( 295 ) for 



the borderline of the support of eigenvalue density yields 
a closed line on the complex plane until pAjj ~ 30, af- 
ter which the line opens from below. This opening is 
reminiscent of the gap predicted by the analytic theory 
for the eigenvalue distribution of the matrix C = ReG 
in Sec. |II.I.4| This signals that an important change in 
behavior might be expected at this density. And indeed, 
we observe that a 'hole' opens in the eigenvalue density 



for pAp > 30. As we see in Fig. 16 c) , this hole is per- 
fectly described by Eqs. (233) and (|234[) which we now 



solve in the biorthogonal basis of right \TZ a ) and left \C a ) 
eigenvectors of the operator T. These eigenvectors obey 
T\lZ a ) = Ha\7tq) and T^\C a ) — ii* a \C a ). In this basis, 



Eqs. (233) and (234 1 read 



9_ 
N 



Si 

a 

1 . , flat 



Ma 



li a u*p(C a \Cp)(1lp\1l a ) 



N tfi 0--Ma)(l-9lifiy 



(301) 
(302) 



where we made use of the fact that TrT = and therefore 
TrS'o = gTrTSo [see Eq. p85] ]. The problem essentially 
reduces to solving the eigenvalue equation 



3 exp(ifc |r-r'|) 
' 1 fcolr-r'l 



TZ a (r') =fi a Tl a (v), (303) 



where r 6 V. As follows from Eq. (288), TZ a (r) is also 
an eigenvector of the Laplacian, A r lZ a (r) = ~K^TZ a (r), 
with n a — K(l//i Q ). In a sphere of radius i?, using the 
decomposition of the kernel of Eq. ( 303 ) in spherical har- 



monics (Gradshteyn and Ryzhik 19801: 



exp(ifc |r - r'|) 
fenlr-r'l 



ji [k mm(r,r')} 



;=o m=-l 

x hf ] [k ma, X (r,r')]Y lm (9^)Y lm (e'^')*, (304) 



it is quite easy to find that (Svidzinsky et al. 2010) 

TZ a {r) = TZi mp (r) = Ai p ji(Ki p r)Y lm (8,(f>), (305) 

where 8 and <f> are the polar and azimuthal angles of the 
vector r, respectively, ji are spherical Bessel functions of 
the first kind, are spherical Hankel functions, Y\ m are 
spherical harmonics, Ai p are normalization coefficients, 
and a = {l,m,p}. Furthermore, coefficients ki v obey 



(Svidzinsky et al. 2010) 



Kip _ ji(Ki p R) h^jkpR) 
k ji-i(Ki p R) h^ikoR) 



(306) 



Integer p labels the different solutions of this equation for 
a given I. Hence, the eigenvalues 



Hp 



1 



2tt 2 (^ p /fco) 2 - 1 



(307) 



are (21 + l)-times degenerate (m g [—1, 1]). 

In the limit kgR — > 00, for / <C kpR and I <C ki p R, we 
can use asymptotic expressions for the spherical functions 



in Eq. (306) to obtain 



ij n / ni p + fc 

2 V Kl p - fc 



-KipR+(~+p)Tt. (308) 



3G 



In this limit, the eigenvalues \xi p are therefore localized 
in the vicinity of a roughly circular line 25 in the complex 
plane given by 



k{1/h) - k 



«(1/m) + k 



= 1. 



(309) 



Let us now study the eigenvectors. Using standard 
properties of spherical harmonics and spherical Bessel 
functions (Morse and Feshbach 1953), we can show that 



(n*t mp \n Vm . p ,) = (-i) m A%^ Mk 1p r) 2 



i? 3 r 

2 



ji-i(ni p R)ji + i(KipR) 



(310) 



From the normalization condition (jO,i mp \TZi' m > P ') = 
Si,l<S m ,m>5p,p', we find that £ imp (r) = (-l) m TZi^ m ) p (r)* 
and 

/~2~ 1 

Alp = \ — (311) 

V R \Jji{KipR) 2 - ji-i{nipR)ji + i(nipR) 
On the other hand, we also have 



(*Rl7np \R\' m' p' 



R 2 Ai p Aip> 
K lp' K lp 



nip>ji-i{Kip>R)ji(K* lp R) 



rfpji-i(n1pR)ji(nip>R) 
Sl,l'S m ,m', (312) 



and {Ci m p\Ci' m 'pi) — (lZi m p\lZi mp >)5ij>5 m:m i . It is now 
convenient to introduce a new coefficient 



C, 



ipp' 



tfpRji- 1 (n*i p R)ji (kip' R) 



Ki p > Rji-i [ni pl R)ji (k* R) 



K 2 p ,R 2 K lp R 2 



-i -l 



jl(KlpR) 2 - jl-l(Ki p R)jl + l(K* p R) 



jl(Klp'R) 2 - j/_l(Kip'i?)j; + l(Ki p /i?) 



in terms of which Eqs. (301 1 and (302) become 

,2 



^^(2J + l) M f p 



1 9 

- + — 

9 i-SA^p 



(21 + l)^ip'^ p C[ P p' 



i p p- 



7 (!- 9lMp>)(l - gvipY 



(313) 



(314) 



(315) 



25 An equation of a circle can be found from Eq. j309| by ex- 
panding in series in 1/pAg. The resulting equation is 
(x + p\l/8n 2 ) 2 + (y-K/2 + l) 2 = K 2 with TZ = 4-y/3W(4:k R), 
fi = x + iy, and W(t) the Lambert function (a function in- 
verse of }(W) = We w ); W(t) ~ In* for |t| > 1. Hence, 
U ~ 47/31n(4fc fi) for k R > 1. 
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FIG. 17 Mean maximum value of the imaginary part of eigen- 
values A of the N x N random Green's matrix G. Analytic 



results (solid lines) following from Eqs. (3141 and (3151 are 
compared with the results of numerical diagonalization for 
three different matrix sizes N (symbols) . Analytic results de- 
pend both on 7 and pAp, except for pAg < 10 when they 
reduce to Eq. (2961 (dot-dashed line). The dashed line rep- 



resents the prediction of the diffusion approximation ( 300 1 . 
The horizontal axis is the on-resonance optical thickness of a 
cloud of cold atoms, relevant for the problem of random lasing 
discussed in Sec. IIV.D.2I 



To find the borderline of the support of eigenvalue den- 
sity of the matrix (282) shown in Figs. 



16 c) and 16 ;d), 



we a pply the following recipe. (1) Find solutions ki p of 
Eq. ( |306 ) numerically and then compute the correspond- 
i ng n i p . (2) Compute the coefficients Ci pp i using Eq. 
(313). ( 3) F ind lines on the complex plane 1/g defined 
by Eq. (315). (4) Transform the lines on the complex 
plane l/g into contours on the complex plane z using 
Eq, 



(314) 



The latter contours are the borderlines of the 
support of eigenvalue density p(A) . 

At high density, the crown formed by the eigenvalues 
blows up in spots centered around fi a — i. wh ere fi a are 
the eigenvalues of T, as we show in Fig. 16 'd). When 
the density is further increased, the clouds of eigenvalues 
of A turn clockwise along the circular line given by Eq. 
(309) and shrink in size. The eigenvalues A eventually 
i. 



become equal to fi a — i. They then fall on the circular 
line (309) and the problem looses its statistical nature. 



As follows from our analysis, the parameter 7 controls the 
overall extent of the support of eigenvalue density V on 
the complex plane, whereas its structure depends also on 
the density pXg. At fixed 7, T> goes through a transition 
from a disk-like to an annulus-like shape, and eventu- 
ally splits into multiple disconnected spots upon increas- 
ing /jAq. The transition from disk-like to the annulus- 
like shape is reminiscent of the disk-annulus transition in 
the eigenvalue distribution of rotationally invariant non- 



Hermitian random matrix ensembles (Feinberg and Zee 



1997a) (see the discussion at the end of Sec. III.C.l). 
Quite remarkably, Eqs. (314) and (315) properly cap- 
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ture the transition to the continuous medium regime 
(high density) and to the small sample regime (small 
koR). To illustrate this point, we calculated (max(ImA)) 
from Eqs. (314) and (315), and found excellent agree- 



ment with numerical results at all valu es o f parameters, 
including high densities pAfj, see Fig. 17 In contrast, 
the analysis of the lower part of the spectrum in Fig. |T6| 
shows that the theory fails to describe it with sufficient 
accuracy, especially at high densities p\q > 10. 



d. Super- and subradiant states. An important additional 
feature of the numerical results in Fig. [16] that is not de- 



scribed by Eqs. (233) and (234) is the eigenvalues that 
concentrate around the two hyperbolic spirals, |A| = 
1 / arg A and its reflection through the origin. These spi- 
rals correspond to the two eigenvalues ±Gi2 of the ma- 
trix (282) for N = 2. The eigenvectors corresponding 



to these eigenvalues are localized on pairs of very close 
points |rj — Yj \ <C A . They are analogous to super- and 



subradiant states of a pair of atoms ( Gross and Harochc 



1982) a nd correspond to the so-called pro ximity reso- 
nances (Heller 1996 Rusek et al. 2000b): the super- 



radiant state corresponds to ReA, ImA > 0, whereas the 
subradiant — to ReA, ImA < 0. Numerical results show 
that in the limit of p\$ — > oo, the lower branch is much 
more populated than the upper one. 

A rough model that partially mimics this behavior is 
given by the N x N matrix: 



G = G 



12 



/0 1 



V 



'•• i 
i oy 



(316) 



where G12 = e 1,Co ' ri ~ r2 '/fco| r i~ r2 1 , and ri and r 2 are ran- 
domly chosen points inside the sphere of radius R. This 
matrix has two different eigenvalues: the non-degenerate 
eigenvalue A = (N — l)Gi 2 corresponds to the superradi- 
ant state (1, . . . , 1) /y/N ; the (N— l)-degenerate eigenval- 
ues A = — G12 correspond to subradiant states localized 
on pairs of points (1,0,..., 0,-1,0, ...,0)/\/2. In the 
limit N — > 00, only subradiant states contribute signif- 
icantly to the eigenvalue distribution of G. Using the 
definition (194) we can show that the latter is then given 
by: 



PK ' (fc i?)3|A|2 ^2fc i?|A| 



S ( argA 



1 

|A| 



(317) 

where s(x) = 1 — 3x/2 + x 3 /2. Loosely speaking, the 
true eigenvalue distribution of the Green's matrix G is a 
superposition of Eqs. ([227|), (f228j) and ([317]) . With the 



qualitative picture of the Dyson gas in mind, we could 
say that the lower 'branch' |A| = —1/ argA plays the role 
of a channel for the gas of eigenvalues, through which the 
latter can escape from the eigenvalue domain predicted 
by Eqs. (233) and (234). This effect is more pronounced 



at high density because the eigenvalues accumulate near 
the line ImA = —1, so that the repulsive interaction be- 
tween eigenvalues forces the latter to flow into the lower 
branch. From numerical results for N < 10 4 , we estimate 
the statistical weight of subradiant states to be impor- 
tant at large densities, of the order of 1 — con st /{pX^Y 
with p ~ 1 (Goetschy and Skipetrov 2011b). This is 



consistent with the estimation of the number of subradi- 
ant states in a large atomic cloud by Ernst (1968). As 
explained earlier, the lack of the spiral branches corre- 
sponding to super- and subradiant states in the theory 
can be traced back to the assumption of statistical inde- 
pendence of the elements of th e m atrix H in the repre- 
sentation A = HTW [see Eq. |Io|]. 

An important implication of the existence of the hyper- 
bolic spiral branches in the eigenvalue distribution of the 
Green's matrix G is that, at least at low density, quanti- 
ties such as (min(ReA)) or (min(ImA)), that are a priori 
difficult to c alculate, can be r eadily found from 2-body 
intera ctions ( Goetschy] |2011[ Skipetrov and Goetschy 
20Tll. 



e. Density of eigenvalues. Equations (231) and (232) al- 



low us to analyze not only the borderline of the eigen- 
value domain, but also the eigenvalue density p(A) inside 
it. Very generally, p(A) is roughly symmetric with re- 
spect to the line ReA — and decays with ImA. In the 
regime of low densities pAg < 1, an approximation of 
Eqs. (231) and (232) can be obtained by replacing the 



operator Si by Sq. This amounts to n eglecting the term 
c 2 TTt in the denominator of Eq. (Eiol. Then, Eqs. (EH) 



and ( 232 ) reduce to two equations in which the reso 



vent 



g(z) and the eigenvector correlator c(z) are decoupled: 



9{z) 



N ±1O 



^TtSqSq 



c(zf = \g(z)f 



N 



TrS Sl 



(318) 
(319) 



Assuming explicitly that the N points are distributed in 
a sphere of radius R, we can make use of the results of 
Sec. |III.D.2.a to compute traces in these equations, so 



that Eqs. (318) and (319) become 



z* -2 1 g(z)*h(iK[g(z)}*R + ik R) 



<zf = \g(z) 



2<yh{2lTnn[g{z)]R) 
1 

~ 2-fh(2lmK[g(z)]R)' 



(321) 



where the functions n{g) and h(x) are defined by 
Eqs. (291) and (294), respectively. We find the resolvent 



g(z) by solving Eq. (320) numerically and then evaluate 
the eigenvalue density p(A) with the help of Eq. (197). 
Note that Eq. (320) applies only within the eigenvalue 
domain T> given by Eq. (295}. Fig. 
distribution p(A) obtained in this way for 
p\g = 1, together with the result of numerical diagonal- 
ization. 



shows the full 
N = 10 4 and 
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FIG. 19 Marginal probability density of the imaginary (left column) and real (right column) parts of eigenvalues A of the 
N x N random Green's matrix ( 282 1 at two different densities p of N points r< randomly chosen inside a sphere of radius 
R. 7 = 9N/S(koR) 2 . Results of numerical diagonalization (blue solid lines) obtained for N 
realizations are compared to the solution of Eq. ( 320 1 (red dashed line) . 



10 after averaging over 10 



f. Marginal distributions of real and imaginary parts of A. 
An interesting link exists between the probability den- 
sities of the eigenvalues of matrices C and S studied in 
Sec. |II.I.4| and |II.I.3[ respectively, and the probability 
densities of real and imaginary parts of the eigenvalues 
of the Green's matrix G. Remember that C and S de- 
fine the real and imaginary parts of G. We note that 
at low densities pAp < 1, the eigenvalue distributions 
of G, ReG and ImG are parameterized by a single pa- 



rameter, their second moment 7 = (|A| 2 ). It is thus 
also the case for the marginal distributions p(ReA G ) and 
p(ImA G ). In addition, Eqs. ( |283[ ) and ( |284[ ) suggest that 
7 = (A| cG ) = 2((ReA G ) 2 ) "(ALg) ^(ImA G ) 2 ), as 
long as koR ^> 1 and the density is not too high. It 
is therefore reasonable to conjecture that p(ImA G ) and 
p(RcA G ) may be described by equations for p(Aj mG ) and 
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p(Ar c g) with 7 replaced by 7/2: 

p(RcA G , 7 ) ~ p(A ReG , 7 /2), (322) 
p(ImA G , 7 ) ~ p(A ImG)7 /2). (323) 

This conjecture was verified by comparing analytic re- 
sults for ReG and ImG with numerical simulations for G 
and turns out to hold with satisfactory a ccuracy as long 
as 7 < 2 (|Skipetrov and Goetschy| |201l|. For 7 > 2, the 



distributions of ImA G and Ai mG were found to be very 
different, whereas the distributions of ReA G and A ReG 
remain similar. 

The marginal probability distributions p(ImA G ) and 
p(RcA G ) can be obtained by projecting p(A G ) following 
from Eq. (320 1 on the real and imaginary axes. As we 
show in Fig. |19[ a good quantitative agreement is found 
with numerical simulations for pAfj < 1. At higher densi- 
ties, Eq. (|320 | is not a good approximation for Eqs. (231 ) 
and ( 232[ ) anymore. At the same time, Eqs. ( |231[ ) and 
(232 1 arc difficult to solve exactly because g(z) and c{z) 
are coupled and TtSiSq has no 'simple' expression in 
the biorthogonal basis of eigenve ctors o f the o perat or T, 
in contrast to TrS Sl [see Eqs. |249| , ( pib] , |314| , and 



is U = -i Yhj /( r * ~ r j)[ A ^( t ) ~ A ^(i)] 2 , wher e the 
second derivative of the pair potential — /(r, — r^) is sup- 
posed to be a rapidly decaying function of its argument. 
Vibrations in this model can be characterized by study- 
ing the Hessian matrix A defined by Eq. ([9| with u = 1. 

In the limit of infinite density p — > 00, this model de- 
scribes an elastic continuum medium. We deduce from 



in Eq. ( 105 ) is given by 



the analysis of Sec. II. F. 2 that the resolvent gk(z) defined 



e(k) 



(324) 



with the dispersion relation e(k) = p[f (k) — /o(0)]. 
In particular, for spherically symmetric smooth /(r) at 
small k we have e(k) = c 2 k 2 . This describes the prop- 
agation of sound with a speed c and a linear dispersion 
relation to = \fe(\s) = ck. Interesting and still not fully 
understood properties arise when the density p is high 
but finite, i.e. when a certain amount of disorder is intro- 
duced in an otherwise homogeneous medium. Attempts 
were undertaken to use the Euclidean RMT to explain 
some of these properties in the framework of the high- 
density expansion developed in Sec. |II.F.2| 



IV. APPLICATIONS 



1. Brillouin peak in the dynamic structure factor 



Let us now see how the mathematical results reviewed 
in the previous sections can be applied to understand real 
physical systems. 



A. Vibrations in topologically disordered systems 

Amorphous solids, glasses and supercooled liquids ex- 
hibit a number of interesting vibrational properties that 
still lack a satisfactory explanation agr eed upon by all 
specialists (Klinger |2010 Philips 1981). Euclidean ran- 
dom 



matrices can be used to model the behavior of 
these topologically disordered systems and thu s to pro 
pose possible explanations of such properties (ICiliberti 



et al. [ 2003; Ga nter and Schirmacher|[2011||Grigera et al. 



2011 



et al 



(2001a 



2001a 2002 2003 2001b I. Following Grigera 



let us model a topologically disordered 
three-dimensional system, — say, an amorphous solid, — 
as an ensemble of N 1 identical particles that har- 
monically oscillate around their equilibrium positions r, 
(2 = 1,..., N). The latter are randomly distributed in a 
large volume V with an average density p = N/V . For 
simplicity, we assume that the displacements A.Xi(t) of 
all particles are collinear and take place along a fixed 
direction that we choose to be the x axis of the refer- 
ence frame. The instantaneous position of the particle 
i is then Rj(t) = r,; + e x Axi(t), where e x is the unit 
vector codirectional with the x axis. The vector nature 
of displacements can be taken into account but does not 
influence the qualitative conclusions of the scalar anal- 
ysis (Ciliberti et al. 2003). The energy of the system 



One of the main quantities measured in inelastic x-ray 
and neutron scattering experiments designed to charac- 
terize the topologically disordered systems is the dynamic 
structure factor (DSF) S(k, u>). For a system of N iden- 
tical particles it is ( Hansen and McDonald 1986 ) 



1 N r 
J M = jyE J 



dte 



iut /pik-IRiW-R^O)]^ 



(325) 



Taking into account only the vibrational modes of the 
system, we can express S(k 7 ui) through the resolvent 
5k (z) defined in Eq. (105) as 



OL-^Tk 2 

S(k,u) = lim Im<7 k (w 2 + ie). (326) 



W7T e->0+ 



Early experiments [see, e.g., (Benassi et al. 1996 



Monaco et al. 1998)] showed that at k < fc , where 
ko is the position of the first maximum of the static 
structure factor, DSF exhibited a low-frequency peak 
with a width T that scaled approximately as T oc k 2 . 
This was very surprising because naive arguments sug- 
gest that Rayleigh scattering in the low-freque ncy (and 



long- wavelength) limit should lead to T oc fc 4 . Grigera 
et aL J 2001a[ ) used the self-consistent equations (120) 
and |mj) to find that Imcr k (z) oc z( d ~ 2 )/ 2 k 2 and T = 
Im<Tk JuF)/uj oc k 2 for d — 3, in contradiction with the ex- 
pectation for Rayleigh scattering but in agreement with 
experimental observations. This result was confirmed by 
a perturbative calcu lation up to the second orde r in Xlp 
(Grigera et al. 2001b). However, more recently Ganter 
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FIG. 20 Reduced density of states v(u)/oj 2 of a topolog- 
ically disordered system at two different densities pa 3 (log- 
log plot). T he s oli d lin es are obt ained by numerical solu- 
tion of Eqs. p20| , (m} and p28| for the ERM defined by 
/(r) = — exp(— r 2 /2a z ) at two densities pa 3 = f and 0.5. 
The dashed lines follow from Eq. < |328| l with g^(z) replaced 
by (?k(z), in the limit of uj — > 0. 



and Schirmacher (2011) pointed out the incompleteness 



of the previous analysis and presented a refined calcu- 
lation of the resolvent gk(z) up to the second order in 
1/p leading to Imerk(z) oc z d l 2 k 2 and hence to T oc k 4 . 



Shortly after that, Grigera et al. ( 2011 ) demonstrated 
the exact cancelation of terms oc z {a ~ ^ 2 k 2 in all or- 
ders in 1/p and also found an additional contribution 
oc z( d ~ 2 ^ 2 k 4 to Imak(z). This again leads to T oc k 4 
and shows that the simple Euclidean RMT model that 
we presented above does not seem to be sufficient to ex- 
plain the anomalous scaling of T with k. Recent experi- 
mental findings (Monaco and Mossa 2009) suggest that 
T oc k 2 is characteristic for the region of relatively high 
frequencies, where strong temperature dependence of T 
is observed and the simple model of harmonically oscil- 
lating particles that we consider in this section does not 
apply. The model becomes justified at lower frequencies, 
where T is independent of temperature, and its predic- 
tion T oc k 4 , which is in agreement with the behavior 
expected for Rayleigh scattering, is verified by the ex- 
periment (Monaco and Mossa 2009[ ). 



2. The boson peak 

Another interesting property of topologically disor- 
dered systems is the so-called 'boson peak' in the density 
of states (DOS) v{u). For our model of a topologically 
disordered system, DOS can be expressed through the 
resolvent (Tl9l) as 



2w .. T , 2 • \ 
lim lmg(w + ie) 

7T £->0+ 



(327) 



An equation for g(z) can be obtained from Eqs. (120) 



and (121) using Eq. (106) ( Grigera et al. 2001a): 



1 



pg{z) 



- + / (0) - Ag{z) 
P 

d3t l J! / ll t \ 

^3/o(q) 2 <?k(*), 



(328) 



where A 
sity (p - 



J d 3 q/o(q) 2 /(27r) 3 - In the limit of high den- 



oo) and low frequencies (Rez 



0), 



trum v(iS) oc u} 2 /p 5 ^ 2 (Grigera et al. 



we can substitute g^j z) in the last term of Eq. (328) by 
,g k (z) defined in Eq. ( 324) and obtain the Debye spec- 

The 



2001a 



2002) 



quadratic scaling of DOS with frequency is precisely the 
behavior that follows from the so-called Debye model for 
DOS in crystals. The disordered nature of an amorphous 
solid manifests itself in an important correction to the 
above scaling at finite p. Experimentally, one observes 
states in excess of the Debye law that give rise to a peak 
in the reduced DOS v(oj) /uj 2 , commonly known as the 
boson peak ( Klinger 2010 1 . The origin of the boson peak 



is still a subject of active research activity that we do not 



intend to review here [see, e.g., ( Chumakov et al. 2011 ) 



for a recent work an d a selection of references] . 
~ (|2001al ~ 



et al. 



Grigera 



2002) argued that the peak in DOS can 



be obtained from the self-consistent equations ( 120 ) and 
(121). Indeed, solving these equations for g^jz ) numer- 
ically, substituting the solution into Eq. (328), solving 



the latter for g(z), and then using Eq. (327), we obtain 
the reduced DOS shown in Fig. 20 At the lower density, 



we observe an excess of low-frequency states with respect 
to the Debye limit shown by dashed lines, even though 
no well-defined peak structure can be identified. Such 
an excess is less important or even absent at the higher 
density. 

The appearance of the boson peak in DOS is be- 
lieved to be generic for a certain class of ERM ensem- 
bles, independent of the specific function /(r) that gen- 
erates the ensemble, as far as the latter decays suffi- 
ciently fast with r. Grigera et al. ([2002) studied it for 
/(r) = (1 — ar 2 /a 2 ) exp(— r z j2a z ), where the parameter 
a were varied. They also checked that computation of 
DOS using the method of moments yields similar results, 
thus proving that the appearance of the boson peak in 
ERM theory is not an artifact of self-consistent equations 
(120) and (121). In general, low- frequency vibrations in 



topological 



y disordered systems are still a subject of in- 



tense theoretical research, as witnessed, for example, by 



the recent works of Amir et al. (2012 2010) 



3. Anderson localization 

A universal phenomenon that is expected to take place 
for waves in disordered systems is Anderson localization 
(Anderson, 1958). In 3D infinite systems, Anderson lo- 



calization manifests itself by a transition from extended 
to localized states when varying either the strength of 
disorder or the energy of the state. 
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Amir et al. ( 


2012 


); 


Ciliberti et al. 


Wu ( 


2009 


);|Krich and Aspuru-Guzik 



(2011 



(2005); Huang and 



studied An- 
derson localization in variants of the ERM model con- 
sidered in this section. Focusing on liquid instantaneous 
normal modes (INMs), |Ciliberti et al. ( |2005[) generated 
random equilibrium configurations of points {rj} in a re- 
alistic system using a previously developed Monte-Carlo 
algorithm at a fixed density and several temperatures. 
These random configurations were used to obtain the 
pair correlation function g(r) that allows to compute the 
probability density function (PDF) of elements of the 
ERM A under assumption that higher-order correlations 
can be neglected. In its turn, this distribution allows 
one to obtain an equation for the PDF of diagonal el- 
ements of a matrix resolvent Q{z) = (z — A) . The 
latter equation was solved numerically using a popula- 
tion dynamics algorithm and the stability of a purely 
real 'population' of Q nn was analyzed. The idea behind 
is that when a small imaginary part is added to Q nn and 
the latter is evolved according to the above-mentioned 
algorithm, lu\Q nn grows in time in the extended phase 
or decays in the localized phase. Indeed, in terms of the 
eigenvectors R n and eigenvalues A„ of the matrix A we 
can write 



z - A 



Im<? nn (z) = Im 



(329) 



where z — A + ie and e — > + . When A belongs to the 
part of the spectrum for which eigenstates are extended, 

\J2^LiK( r i) R m( r i)\ 2 ~ 1 / N and Im £™ ~ 1- In con- 
trast, for A in the part of the spectrum corresponding to 

localized states, i?* (rj)i? m (rj) will be sizable only for a 
small fraction of sites and ImQ nn = for most n. This 
approach allowed Ciliberti et al. ( 2005[ ) to show that two 
mobility edges exist in the spectrum of the matrix A (one 
at negative and another one at positive A) and to roughly 
determine their positions. 

A different approach to essentially the same prob- 
lem was employed by Huang and Wu (2009). They 



were able to differentiate between the parts of the eigen- 
value spectrum corresponding to extended and local- 
ized states by analyzing the level spacing statistics, i.e. 
the statistical distribution P(s) of normalized spacing 
s = (A„ +1 — A ra )/(A n+1 — A„) between ordered eigen- 
values A„. This distribution is expected to have different 
shapes in the extended part of the spectrum, where re- 
pulsion of eigenvalues is expected [P(s) cx s for small a], 
in its localized part, where P(s) cx exp(— s), and at the 
mobility edge, where P(s) takes an intermediate form. 
Using extensive numerical simulations, |Huang and Wu 
(2009) not only located the two mobility edges in the 



spectrum of ERM A quite precisely, but was also able to 
determine the critical exponent of the localization tran- 
sition from the finite-size scaling of the second moment 
of P(s). The value of the critical exponent v ~ 1.6 is in 
agreement with the one found for different models of the 
same universality class in 3D. 



A similar estimate of v was obtained by Krich and| 
Aspuru-Guzik (20111 who studied ERMs generated by 



an exponentially decaying function /(r) (see Sec. II. 1. 2) 
They performed a finite-size scaling analysis using power- 
ful numerical methods. An analytic treatment of Ander- 
son localization in this model was developed very recently 
by Amir et al. ( 2012[ ) who found the critical frequency of 
the localization transition as a function of density of ran- 
domly distributed points and the dimensionality d of 
the system. 



B. Electron glass dynamics 

An electron glass is a highly disordered solid in which 
most electronic states are localized and where long-range 
Coulomb interactions play an important role. It is re- 
ferred to as a glass because it exhibits glassy behaviors, 
such as slow relaxation of the conductance and depen- 
dence of the relaxation on the time of application of a 



perturbation. The latter property is called 'aging' (Amir 



et al. 2011) 



A simple model of the electron glass is an ensemble 
of exponentially localized electronic states randomly dis- 
tributed in space and weakly coupled by phonons. At 
equilibrium, the hopping rate between the states local- 
ized at Yi and Yj is given by the Fermi's golden rule 

7° cx f[Ei)[l - fiE^e-^il + b(Ei - Ej)}, (330) 

where / and b are the Fermi-Dirac and Bose-Einstein dis- 
tributions that characterize the statistics of electrons and 
phonons, respectively. E { = e l + J2j 7 ti e2 f( E j)/ r ij > E j 
is the potential energy of the state i and ej is the en- 
ergy in the absence of Coulomb interactions. If Ei < Ej, 
the expression in square brackets in Eq. ( 330 ) have to be 
replaced by b(Ei — Ej). Within the local mean- field ap- 
proach, the average occupation number nt evolves, under 
a small perturbation, as drii/dt = Y^jilji ~ 7y); where 
7i ? - is obtaine d from Eq. ( |330[ ) upon replacing f(Ei) by m 
(Amir et al. 2008). In the vicinity of a metastable state 



characterized by the electronic configuration n , the lin- 
earized equation of motion for <5n = n — no takes the 
form of a master equation d<5n/di = ASn, where the off- 
diagonal elements of the matrix A are 



A, 



n§(l-n§) 



kjtj,i 



1 



1 



(331) 



J2ijii^ij guarantee 



and the diagonal elements An — 
the particle number conservation. 

Near a stable point, the eigenvalues of A are negative 
and their statistical distribution p(A) determines the re- 
laxation p roperties of t he gla ss. A numerical study per- 
formed by Amir et al. ( 2008 ) shows that for parameters 



relevant to experiments, p(A) is similar to the distribu- 
tion obtained by keeping only the first term in Eq. (331 ) 



and neglecting the remaining energy dependence, such 
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that Aij is replaced by Am — e Therefore, the ex- C. Waves in open chaotic systems 



ponential ERM A studied in Sec. |II.I.2| appears to be of 
fundamental importance for the study of slow relaxation 
processes in an electron glass. In the regime of low den- 
sity p(, d <C 1, where almost all eigenvectors of A are local- 
ized, one readily finds from Eq. (178) that p(A) cx 1/|A|. 
This behavior gives rise to a logarithmic relaxation of a 
small deviation 6n(t) from the metastable state. Indeed, 
expanding 5n(t) over the basis {Rfc} composed of eigen- 
vectors Rfc of the matrix A: Sn(t) — J2k c feR'fe e_ ' Afc ' t ) 
and assuming that the eigenvectors are uniformly excited, 
we get 



\Sn(t)\ 



-Al 



dA- 



A 



-7b - ln(A min i) (332) 



for 1/A max < t < 1/A m j n , where je is the Euler constant 



(Amir et al. 2008). It is worth noting that the Coulomb 



interaction does not play any role in this slow relaxation, 
since the latter has been neglected in the replacement of 
A by A. 

In real experiments, a logarithmic decay is observed 
not directly for the deviation \Sn\ of the occupation 
number but for the conductance 5a of the electron-glass 
(Vaknin et al. 2000). It is believed that the relation 



<5cr(|<5n|) is linear (Amir et al. 2008[ ), which implies that 
a kick out of equilibrium increases the conductance before 
the latter starts relaxing slowly. One possible explana- 
tion for this effect is related to the long-range Coulomb 
interaction as follows. At low temperature and at equilib- 
rium, the interaction creates a soft gap in the density of 
states at Fermi energy, known as the Coulomb gap, that 
reduces the conductance. When a small perturbation is 
applied, the Coulomb gap is suppressed, causing the con- 
ductance to increase before relaxing back to equilibrium 
(Amii etal. 2011). Hence, Sa(t) evolves according to Eq. 



(332). This result allows us to understand aging experi- 



ments, in which a gate voltage is applied during a time 
t w to a sample initially at equilibrium. By taking into 
account the fact that not all the modes have the time to 
relax during t w , it is easy to show that at time t (mea- 
sured from the end of the excitation), the conductance 
relaxes as 



5a(t,t w ) =C\n{l + t w /t), 



(333) 



where C is a nonuniversal constant (Amir et al. |2009 1. 
The fact that Sa(t, t w ) depends on the ratio t w /t is known 
as 'full aging' and is a direct co nsequ ence ofp(A) ~ 1/|A|. 
The behavior predicted by Eq. ( 333 ) was observ ed in elec- 
tron glasses su ch as InO (Vakni n et al. 2000) or granu- 
lar aluminium (Grenet et al. 2007), as well as in struc 



tural glasses such as the plastic Mylar (Ludwig et al. 
20031. This suggests that the ERM model of relax- 



ation described in the present section might be a general 
paradigm for aging in glasses of different physical origins 
( |Amir et al. [|2009"] ). 



In classical mechanics, a chaotic system is a determinis- 
tic, often quite simple system that exhibits extreme sensi- 
tivity to initial conditions. This means that in practice, it 
is impossible to make any predictions for the state of such 
a system after a sufficiently long time. Billiards having 
irregular shapes are examples of classically chaotic sys- 
tems. Interestingly enough, when a quantum (or wave) 
problem is solved in a system that is classically chaotic, 
the chaotic behavior is suppressed and the evolution of 
the wave function of a quantum particle (or of the am- 
plitude of a classical wave) can be calculated without 
problems. However, the resulting solution bears subtle 
signatures of the chaotic behavior that the system would 
exhibit if it were classical. These signatures are subject 
of the field of 'quantum chaos' (Haake 2010). 



Particularly interesting physics takes place in open 
chaotic systems that couple to the exterior world through 
M 'channels'. The statistical properties of the matrix 
X considered in Sec. |III.D.1| are strongly reminiscent 
of those of effective Hamiltonians used to characterize 
open chaotic systems (|Fyodorov et al.\ 2005 |Fyodorov 
and Sommers[ |1997| |Mahaux and W cidcnm uller| |1969| 
Mitchell et al. . 2010 1. If we recall the physical meaning of 



the matrix X, we understand that this analogy is not an 
accident. Mahaux and Weidenmiiller ( |1969[ ) introduced 
the random matrix model for the M x M scattering ma- 



trix of an open chaotic system [see also (Mitchell et al. 
2010])]: 



S(E)=I M -\arf 1 H, 

& — HcS 



H, 



eff 



—HH\ 
2 



(334) 
(335) 



where Hq is a Hermitian matrix that describes the closed 
part of the system under consideration, E is the energy 
of the incoming wave, H is a N x M matrix that con- 
tains entries coupling the N internal states to the M 
external channels, and a > is an overall coupling pa- 
rameter controlling the 'degree of non-Hermiticity' of the 
effective Hamiltonian _ff c ff- Hq is commonly drawn from 
the Gaussian ensemble Q, and H is chosen such that 
is a Wishart matrix. Randomness in Hq and HH^ 
is assumed to be independent, meaning that Hq and 
HW are asymptotically free matrices. The eigenvalue 
distribution of -ff e ff was considered previously by |Haake 
et al. ( 1992 ) (with the help of the replica trick) , [Lehmann 
et al. ( 1995 ) (using the supersymmetry method) and 
Janik et al. (1997a]) (using the free probability theory). 
The splitting of the domain of existence of eigenvalues 
in two parts was observed when a was increased. This 
is slightly different from the matrix X studied in Sec. 
III.D.l that has elements with equal variances j/N of 
real and imaginary parts (hence always the same degree 
of non-Hermiticity). 

To fully understand the similarities and the differences 
between the two models, let us consider a realistic MxM 
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scattering matrix S(lu) describing the scattering of light 
by N randomly placed atoms with polarizability a(uj) 
(|Goetschyl[20Tll): 



S(u) = hi + TH^ 



1 



Aiuj)- 1 - G 



H, 



G 



HTH 



(336) 
(337) 



where uj is the frequency of light, A(uj) — 
diag[/cQa(u;)/47r] is the polarizability matrix, H is defined 
by Eq. pTj) , and G is the Green's matrix considered in 
Sec. |IILD.2| By choosing a simple model for the atomic 
polarizability a(u>) = -(47r/fc^)(r /2)/(w - w + iT /2), 
we obtain 

S(u) = Iu-Ttf To/ £ H, (338) 

W — il c ff 

h cS = uj i n - y RcG ~ IT [ImG + In] ■ (339) 

The main difference between the effective Hamiltonian 



( 339 1 with respect to Eq. ( 335 1 is that its Hermitian and 



anti-Hcrmitian parts are correlated because they both de- 
pend on the matrix G. A model with uncorrelated Hermi- 
tian and anti-Hermitian parts can be obtained by replac- 
ing G by the matrix X studied in Sec. |III.D.f | However, 
we see from the comparison of the eigenvalue densities 
of matrices X (Fig. [f5| and G (Fig. 16 1 that such a re- 
placement can be justified only at very low densities and 
7 <C 1, when eigenvalues of both X and G are uniformly 
distributed within a circle on the complex plane. 26 If wc 
restrict our consideration to low densities and 7 < 1 and 
accept to replace G by X in the effective Hamiltonian 
(3391, the latter becomes very close to Eq. (335): ReX 



is well approximated by a Gaussian matrix, mimicking 
therefore the random Hamiltonian Hq of Eq. (3351, and 
the anti-Hermitian part of Eq. ( 339 ) S — ImX + In is 
well approximated by the Wishart matrix jHH^ , which 
coincides with the anti-Hermitian part of Eq. ( |335 ) with 
a = 7 = N/M. The eigenvalue domain of X [and hence 
of i/eff of Eq. ( 339 )] splits in two parts upon increasing 
7 (see dashed contours in Fig. 15 1, exactly in the same 
way a s does the eigenvalue domain of H e s defined by Eq. 
fl335| (|Haake et al] |1992[ |Janik et al] [i~997a| |Lehmann 



et al] [1995 1. However, it is important to keep in mind 
that this similarity stems from the approximations made 
for the matrix X that are not justified at large 7 at which 
the splitting takes place. Indeed, we see from Fig. [15] 
that the numerically computed eigenvalues do not split 
in two groups, in contradiction with the prediction of 
the approximate theory shown by the dashed line. The 
splitting of the eigenvalue domain is therefore an arte- 
fact of approximations. The similar splitting discussed 



26 Note that even at low densities and 7 < 1, the eigenvalue dis- 
tributions of matrices X and G exhibit differences due to the 
eigenvalues that fall outside this circle: they follow two hyper- 
bolic spirals for G but form a 'cross' for X. 



by Haake et al. ( 1992 1; Janik et al. ( |1997a ); Lehmann 
et k|(|1995[);|Mitchell et a/.|(|2010|) might also be an arte- 



fact of the model ( 335 ) that fails to describe the correct 
effective Hamiltonian of the system at large a. 27 



D. Waves in random media 



The effective Hamiltonian ( 339 1 for a scalar wave in an 
ensemble of randomly distributed, resonant point scatter- 
ers (atoms) can be rewritten as 



H c s — [ ujq — i 



.Tr 



N 



-G. 



(340) 



We thus see that the Green's matrix G defined by Eq. 
( 282 1 essentially plays the role of the effective Hamilto- 



nian for light scattering in an ensemble of point scatter- 
ers. An eigenvector of G associated with an eigenvalue 
Afc corresponds to a quasi-mode of the system that has an 
eigenfrequency Uk — loq — (Fo/2)ReAfc and decays with a 
rate r fe /2 = (r /2)(l + ImA fe ). Hence, the study of the 
statistical properties of A^ allows us to better understand 
scattering of waves among a large ensemble of scattering 
centers. In this context, a few particular problems will 
be discussed below. 



1. Collective spontaneous emission of large atomic ensembles 

An interesting problem of modern quantum optics con- 
cerns relaxation of elementary excitation in large atomic 
systems. To put it simple, a photon is stored in an en- 
semble of (cold, meaning immobile) atoms, and one is 
interested in the properties (frequency, direction of prop- 
agation, etc.) of the photon re-emitted by the atoms at 
a later time. It is a specific case of the superradiancc 
protocol, with only one photon and no restriction con- 
cerning the size of the system. Theoretically, this prob- 



lem was addressed a long time ago by Ernst (1968), but 



has been popularized only very recently by the group 



of Scully (|Das et al] 2008 


Scully 


2009; 


Scully et al. 


2006[ Scully and Svidzinsky[ 


2009, 2010; Svidzinsky and 


Chang[ 2008 |Svidzinsky et al.\ 2008 |2010| Svidzinsky 


and S cully [ 


2009 


), as well as by Manassah and Fried- 



therein] . The reason for this renewed interest is probably 
the recent development of experimental setups free from 
effects that may obscure cooperative phenomena (e.g., 



27 Note that splittings of the eigenvalue domains are, however, ob- 
served for ERMs S (Sec. [ILL3J, C (Sec. [ILL4} and G (Sec 



|III.D,2 1 . The crucial point is that they do not occur at 7 
N/ (koL) 2 ~ 1, but at koL ~ 1. In particular, in the limit of 
very small sample k^L 1, the cloud of eigenvalues of G with 
the large st imaginary part describe s the phenomenon of super- 
radiance (Gross and Haroche 1982k 
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Doppler effect or near field atom-atom interactions), us- 
ing either cold atoms or ultrathin solid-state samples 
(IRohlsberger et aZ.l 120101). 



Consider N immobile two-level atoms at random posi- 
tions in a three-dimensional space. Each atom has 
the ground state \g) and the excited state |e). The atoms 
can interact with a single quasi-resonant photon of fre- 
quency lo close to the frequency ujq of the atomic tran- 
sition. The quantum state of the whole system can be 



writt en as ( Svidzinsky et al. 2010 Svidzinsky and Scully 
20091) 



N 



!*(*)} = ^2mW-l):g,j:e)\0) 



]T 7k (;)|iV:g>|k) 



k 

JV 

+ :: V :'/.' :<../ :< k . (341) 

i<j k 

where the first sum is over states with one atom (atom 
j) in the excited state and no photons, the second one is 
over states with all atoms in the ground state and a pho- 
ton in the mode |k), and the last sum is over states with 
two atoms (i and j) excited and a photon in the mode 
|k). The last term — in which the photon is virtual — is 
necessary to describe non-resonant processes and recover 
the proper dynamics of the system as a whole. The evo- 
lution equation for the vector f3 = (/3i, . . . , /3/y) reads 



(Svidzinsky et al. 2010 Svidzinsky and Scully 2009) 



d/3 



= -/3(<) + iG/3(t), 



(342) 



where the time is in units of the lifetime of t he excited 
state Tq 1 and G is the Green's matrix (282 1. Accord- 



ing to this equation, and as we alre ady obtained above 
from the effective Hamiltonian ( 340 ) , the real and imag- 
inary parts of an eigenvalue of G yield the frequency 
shift and the modification of the decay rate of the cor- 
responding eigenstate. The frequency shift originates 
from non-resonant contributions (so-called 'off-shell pro- 
cesses') in the light-matter interaction, and is sometimes 



called 'collective Lamb shift' ([Scully 2009; 


Scully and 


Svidzinsky[ |2010||Svidzinsky et al.\ 2010||Svid 


zinsky and 


Scully 20091. Both the decay rate and the 
shift were studied by | Svidzinsky and Scully | 


frequency 
2009) and 



Svidzinsky et al. ( 2010[ ) in the limit of a very dense atomic 
cloud (/oAq ~ oo) by re placin g summation by integration 



in the last term of Eq. (342 1, [G(3(t)] 



This is equivalent to averaging Eq. ( 342 1 over all possible 



configurations {r^} of atoms and amounts to the neglect 
of the statistical nature of the initial problem in which 
both the decay rates and frequency shifts were random 
quantities, dep ending on the positions {r^} of a toms. As 
a consequence, Svidzinsky and Scully (2009) and Svidzin- 
sky et al. (20101 find deterministic eigenvalues A fc . Be 
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FIG. 21 Distribution of dimensionless decay rates y = 
ImA + 1 of eigenstates of an ensemble of TV = 10 4 identi- 
cal two-level atoms in a cube of side L. Results at densities 
pAg = 1, 10, 20, 40, 60 and 100 (curves from top to bottom) are 
compared with the asymptotic law 1/y shown by the dashed 
line. [Reproduced from (Skipetrov and Goetschy 20111.] 



i mpor tance of which was already pointed out by [Ernst 
([19681, are lost. 



In this context, the study of the random Green's ma- 
trix reviewed in Sec. |III.D.2| appears to be quite useful. 
In fact, it allows us to immediately obtain analytical re- 
sults for the statistical distributions of the dimensionless 
decay rate y = 1 + ImA and frequency shift x = ReA 
at low density and small on-resonance optical thickness 
6o = I6 7/3 ( Skipetrov and Goetschy 2011[ ). According 
to Eqs. (322) and (323), p(y) is given by the Marchenko- 



Pastur law ( 181 ) with A replaced by y, while p(x) follows 



from the more involved Eq. ( 187) after replacing A by x. 



In both cases, 7 should be replaced by 7/2 and 60 is the 
only parameter of the distribution. At higher densities 
or f or b p > 1, the more advanced results presented in 
Fig. 



19 



apply. For /jAq > 10, no analytical results for the 
distributions of x and y are available. Numerical simu- 
lations show that the distribution of dimensionless decay 
rates y = ImA + 1 tends to 1 /y as the density is increased 



(Skipetrov and Goetschy |2011 ), see Fig. 21 



In a recent work, Akkermans et al. ( 2008 1 claimed that 



statistical properties of decay rates y can be understood, 
at least qualitatively, by dropping the real part of the 
Green's matrix, inasmuch as the latter is expected to be 
responsible for the collective Lamb shift. This picture 
might not be entirely correct, because the decay rates 
and frequency shifts are related to the imaginary and 
real parts of the eigenvalues of the matrix G, and not to 
the imaginary and real parts of the matrix itself. Nev- 



sides, all subradiant states of the Green's matrix, the 



ertheless, the results of Sec. |III.D.2.f| suggest that under 
certain conditions, the distributions of imaginary parts of 
eigenvalues of matrix G and of eigenvalues of ImG have 
indeed some similarities. 
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y levels of 


a two-level atom in the field of 
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a quasi-resonant pump beam (amplitude Q p , frequency uu p ). 
(b) Energy levels of a three-level atom pumped incoherently 
through the auxiliary level a; F ae 2> To, Foilp. 



2. Random laser 

The emission of light by a conventional laser requires 
two important ingredients: an active medium that am- 
plifies light and a feedback mechanism that ensures that 
light passes through the active medium many times ( Sieg- 



man 



1986 ) . The feedback is typically provided by an op- 
tical cavity that also plays a major role in the selection 
of modes in which the laser emits. In a random laser, 
the cavity is absent and the feedback is ensured by the 



multiple scattering of light (Cao 2005 Wiersma| 2008 1 . 



Euclidean RMT was used to describe the simplest ran- 
dom laser — a random ensemble of atoms under an ex- 
ternal pump (Goetschy and Skipetrov, 2011a). Con- 



sider a gas of N three-level atoms at random positions 
Yi (i — 1,...,N) in free three-dimensional space. The 
atoms are illuminated by a strong external pump field 
resonant with the transition from the ground state \gi) 
to the upper auxiliary level |ai) of each atom. The atoms 
rapidly decay to the upper level |ej) of the laser transition 
at a rate r ae 3> r eg = Tq r ag , as we illustrate in Fig. 
[22|b). Interaction of atoms with the electromagnetic field 
which is near-resonant with the transition from |ej) to \gi) 
(energy difference Hluq) is described by 5iV equations of 
motion for atomic operators that are coupled to the quan- 



tum propagation equation for the electric field ( Goetschy 



2011). After elimination of the electric field, and in the 
semiclassical approximation, these equations can be re- 
duced to a system of 2N equations for the expectation 
values and IT of atomic raising and population imbal- 
ance operators = leMgA and IT = |ej)(e t | - \gi){gi\, 



respectively (Goetschy and Skipetrov 2011a 



d^ 
dt 



dlT 
dt 



i o i 



= -(l + WJILi + Wi-l 



(343) 



21m 



&t Sj 



(344) 
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FIG. 23 Fhe domain T> a (hatched) spanned by 1/a and the 
domain 2?a (blue area delimited by the solid line) occupied 
by the eigenvalues A of the random Green's matrix ( 282 1, (a) 



Incoherent gain, corresponding to the atomic level structure 
of Fig. 22 b). (b) Coherent Mollow gain, corresponding to 
the atomic level structure of Fig. 22 a) with A p — (u p — 
^o)/ro = 1. Lasing starts when T> a and T>a overlap: regions 
(1), (2), (3). The borderline of V A is given by Eq. ([347} with 
the optical thickness bo = 40 in (a) and bo — 140 in (b). 
The dashed lines show the borderline of Da following from 
the diffusion approximati on [Eq . ( 300 )] . [Reproduced from 
( Goetschy and Skipetrov 2011al.] 



rate that determines the equilibrium population imbal- 
ance II° q = (Wi - l)/(Wi + 1). We also recognize the 
random Green's matrix G that couples different atoms. 

The lasing threshold can be found by analyzing the 
linear stability of the stationary solution IT = n° q , Sf — 
of Eqs. (343) and (344). For a uniform pump Wi 



lasing starts whenever 
2W 



l + W 



ImA„ > (1 + W) + ImA, 



W, 



(345) 



where A„ is an eigenvalue of G. This condition can be 
rewritten in a more general form, applicable to any point- 
like scatterers, as 



A, 



1 



a{u) L y 



(346) 



where a(uj) — (fcQ/47r)a(a;) and a(ui) is the polarizabil- 
ity of the scatterer (atom). The lasing frequency lul is 
to be determined from this equation. In Fig. 23 we illus- 
trate the solution of this equation for a three-level atom 
that we considered above [Fig. [23ja)], but also for a two- 
level atom under a quasi-resonant pump [Fig. 23 'a) , see 
the level structure in Fig. 22 ^a)]. In the latter case, the 
polarizability a{uS) is given by the famous Mollow result 
(|Mollow[[T972|. According to Eq. ([346}, to find the lasing 



Here the time t is in units of r o and Wi is the pumping sity pX 



threshold, one has to find an overlap between the eigen- 
value domain T>\ of the Green's matrix and the domain 
T> a spanned by l/a(w) when the parameters of a(io) are 
varied. The two domains 2?a and T> a are shown in Fig. 
23] in blue and gray hatched, respectively, for the param- 
eters that correspond to the random laser just above the 
threshold in both cases. 

As we discussed in Sec. |III.D.2| at a moderate den- 
10, the eigenvalue domain T>\ consists of two 



4G 



parts: a (roughly circular) 'bulk' and a pair of spiral 
branches. Depending on the particular model of atomic 
polarizability a, either the bulk or the branches may 
touch T> a and give rise to the laser effect. The lasing 
threshold due to the bulk of eigenvalues is readily ob- 
tained by combining the analytic equation ( 296 ) f or th e 
borderline of Pa at low density pAg < 10 and Eq. (346): 



^b \a\ 2 h (-b Ima j = 1, 



(347) 



where h(x) is given by Eq. (294). Note that the thresh- 
old depends only on the on-resonance optical thickness 
bo but not on the density of atoms pAg. Interestingly 
enough, for both gain mechanisms represented in Fig. 
23 the threshold condition ( 347 1 involves the eigenvalue 



with the largest imaginary part. The analytic calculation 
of (max(ImA)) presented in Sec. III.D.2.C is in excellent 
agreement with numerical results, see Fig. 17 It is quite 
remarkable that the agreement is present at all values of 
parameters, including high densities pAg 3> 1 that were 
necessary to reach large optical thicknesses frg 3> 1 in 
numerical calculations with moderate N < 10 4 . Because 
it is (max(ImA)) that controls the laser threshold, we 
conclude that ERM theory applies to random lasing in 
ensembles of atoms all the way from weak (pAg <C 1) to 
strong (pAg ^> 1) scattering regime. 



Goetschy and Skipetrov (2011a) and Goetschy (2011) 
have shown that Eq. ( |347 1 predicts a lower threshold for 



lasing than the diffusion approximation, based on Eq. 
(300). They also analyzed lasing due to the eigenvalues 



that belong to spiral branches shown by solid black lines 
in Fig. [23] and that may enter into play in the case of 
Mollow gain, corresponding to the level structure of Fig. 



22 'a) . ERM theory of random lasing seems to be partic- 
ularly adapted to the description of lasing in cold atomic 
systems and may be the method of choice for the descrip- 



tion of recent experiments (Baudouin et al. 2013). 



3. Anderson localization in open random media 



The methods described in Sec. IIV.A.3I and aimed at 
the study of Anderson localization in an infinite medium 
cannot be applied in open random media, where exci- 
tations can leak outside the medium. This is the case, 
for example, for light or sound in random ensembles of 
small scatterers. As follows from Eq. (340), for point-like 



scatterers the study of Anderson localization reduces to 
the study of statistical properties of eigenvectors of the 
random Green's matrix G. According to the Ioffe-Regel 
criterium of localization, in a 3D medium the localization 
transition is expected for ki ~ 1. Here k is the wavenum- 
ber of the wave and I is the scattering mean free path due 
to disorder. None of these is actually known when the 
scattering is strong, but we can estimate ki by replacing 
k by fcg and I by £o = k^/Awp, the on-resonance mean 
free path in the independent scattering approximation 
(Lagendijk and van Tiggelen] 1996). This leads to the 



conclusion that in a system of point scatterers, Anderson 
localization can be expected at pAg > 20. 



a. Distribution of decay rates. First attempts to use ran 
dom Green's matrices to study Anderson localization 
were undertaken by Rusek et al.\ (2000a) and Pinheiro 
et al. (2004). In particular, in the latter work the au- 



thors looked for signatures of Anderson localization in 
the probability distribution p(y) of dimensionless decay 
rates y = ImA + 1 of eigenstates that we defined in Sec. 



IV.D.l Based on heuristic arguments and extensive nu- 
merical simulations, they concluded that p{y) takes a 
universal form p{y) oc l/y (for a certain intermediate 
range of y's) when the regime of Anderson localization 
is reached. Although Fig. 21 seems to confirm this con- 
clusion, one may note from this figure that p(y) tends to 
decay as l/y at densities as low as pAg = 1, which is far 
from the mobility edge expected at pAg « 20. To get a 
deeper insight into this issue, we show cuts of p(A) along 
the imaginary axis ReA = in Fig. [24|a). We clearly 
observe that p(ReA = 0,1mA) decays as l/(ImA + 1), 
even though the density of points pAg is too low to bring 
the system to the localization transition. For 7 < 1, al- 
though p(A) oc l/(ImA + 1), the marginal distribution 
p(ImA) follows the Marchenko-Pastur law [see Fig. 19 
due to the circular shape of the support of p(A). Inci 



dentally, we see that the Marchenko-Pastur law (181) for 



7 < 1 can be seen as a projection of a two-dimensional 
distribution p(A) on the imaginary axis, provided that 
p(A) oc l/(ImA + l) inside a circle of radius y/2^y centered 
at (0, 7/2) and p(A) — elsewhere. The power-law decay 
becomes visible in the marginal distribution p(ImA) (see 
Fig. 21 ) only when the support of p(ImA) is sufficiently 
wide, i.e. for 7 > 1. Because the condition 7 > 1 can 
be obeyed at any, even very low density by just increas- 
ing the number of scatterers N, it seems that no direct 
link can be established between the power-law decay of 
p(ImA) and Anderson localization. This also seems to be 
confirmed by the calculation using the analytic Eq. ( |320[ ) 
at large values of N that are inaccessible f or n umerical 
simulations and low density pAg 



0.1 [Fig. 24 "b) 



b. Properties of eigenvectors. A more direct way of study- 
ing Anderson localization is to look at the eigenvec- 
tors of the Green's matrix. At high enough density 
pAg, the eigenvectors of t hree d ifferent types coexist 



(Goetschy and Skipetrov 2011a), as we illustrate in 



Fig. 25 The eigenvectors of type (d) appear only at 
pAg > 20 and are localized due to disorder. A quantita- 
tive measure of degree of localization of an eigenvector 
R n = {i?„(i"i), . . . , R n (rjy)} can be obtained by comput- 
ing its inverse participation ratio (IPR) 



IPR ri 



(348) 
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FIG. 24 (a) Cuts of the eigenvalue density p(A) of the TV x TV Green's matrix ( |282[ | along the imaginary axis ReA = 0. TV = 10 4 
points Ti are randomly chosen inside a sphere of radius R; 7 = 9TV/8(fco-R) 2 . Results of numerical diagonalization (symbols) are 
compar ed w ith the solution of E q. j3 2Cjn (solid red lines), (b) Marginal probability density of the imaginary part of eigenvalues 



of Eq. (2821. Solutions of Eq. (320i(solid lines) at TV = 10 4 (7 = 0.34), 10 6 (7 = 1.6), and 10 s (7 = 7.4) for pAg = 0.1 are 



compared with the asymptotic law l/(ImA + 1) (dashed line). 




FIG. 25 (a) Eigenvalues A of a single random realization of 
the Green's matrix G (dots) for a cloud of optical thickness 
60 = 40, composed of TV = 10 3 atoms, (b)-(d) Intensities 
\R\\ 2 corresponding to the mode in the subradiant branch, 
localized on a pair of atoms (b), the mode with the largest 
ImA (c) and the mode corresponding to the smallest |A| (d). 



A mode R n = {Rn, R n 
centered at positions of atoms r; 
lx (b), lOOx (c), and 10 x \Rt 



Goetschy and Skipetrov (2011a). 



, R n } is represented by spheres 
and having radii equal to 
2 (d). [Reproduced from 



n 






-15 -10 



10 15 



FIG. 26 (a)-(c) Density plots of the logarithm of the aver- 
age inverse participation ratio of eigenvectors of the Green's 
matrix ( 282 1. For each of these plots, we found eigenvalues of 
10 different random realizations of 10 4 x 10 4 Green's matrix 
numerically (with points randomly chosen inside a sphere 
of radius R), computed their IPRs using Eq. (348 1, and then 
determined IPR(A) by integrating Eq. ( |349| over a small area 
(AA) 2 around A, for a grid of A's on the complex plane, (d) 
Density plot of the logarithm of the eigenvalue density of the 
Green's matrix (2821. The solid red line represents the bor- 



derline of the support of eigenvalue density following from 
Eqs. ((3141 and ((3151. 



An eigenvector extended over M points is characterized 
by IPR ~ 1/M. The average value of IPR corresponding 
to eigenvectors with eigenvalues in the vicinity of A can 
be defined as 



IPR(A) = ^ X y(|]lPR n( 5 2 (A-A n )), (349) 



\n=l 



where averaging is over all possible configurations of TV 
points in a sphere. The numerical analysis of the average 



IPR defined by this equation reveals the following sce- 
nario ( (Goetschy and Skipetrov 2011b). At low density 
pXq < 10, IPR ~ 2 /TV for all eigenvectors except those 
corresponding to the eigenvalues that belong to spiral 
branches [see Fig. 16 'a) and (b) and section |III.D.2.d 



for which IPR ~ |. These states are localized on pairs 
of points that are very close together and correspond to 
proximity resonances (Rusek et al. 2000a). The prefac- 
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tor 2 in the result for IPR of extended eigenvectors is 
due to the Gaussian statistics of eigenvectors at low den- 
sities. For pXl > 10 [Fig. fffl a) and (b)], IPR starts to 



grow in the vicinity of A = a nd reaches maximum val- 
30 [Fig . 26 c)]. Contrary to common 



0.1 at p\% 



ues 

belief ( |Rusek et aL[|2000a| ), neither localized states nec- 
essarily have ImA close to —1, nor states with ImA ~ — 1 
are always localized, as can be seen from Fig. 26 'c). For 
pAp > 30, the localized states start to disappear and a 
hole opens in the eigenvalue density. As can be seen from 
the comparison of Fig. 26 'c) and (d) , it is quite remark- 
able that the opening of the hole in p(A) [Fig. 26 ]d)] pro- 
ceeds by disappearance of localized states [i.e. of states 
with IPR > l/N in Fig. [26^)]. Although this suggests 
a link between the hole in p(A) and Anderson localiza- 
tion, further work is needed to arrive at any definitive 



conclusions (Goetschy and Skipetrov 2011b) 



V. CONCLUSION AND PERSPECTIVES 

In this paper, we tried to give an overview of the cur- 
rent state of the art of the Euclidean random matrix 
theory. We presented a number of approaches that can 
be used to deal with ERMs and gave examples of the 
practical use of ERM theory to understand real physical 
systems. It is clear, however, that much remains to be 
done in this field. Our understanding of mathematical 
properties of ERMs is far from being complete and their 
applications could certainly be much wider than the few 
examples that we considered. 

One possible direction of research in which progress 
is needed concerns the representation of an arbitrary Eu- 
clidean matrix as A = HTH^ that we heavily exploited in 
this review. The assumption that the matrix H has i.i.d. 
elements Hi a allowed us to progress in the calculation but 
it is certainly not sufficient. Correlations between Hi a 
could be taken into account by anal ogy with the work on 
correlated Wishart matrices HW (Marchenko and Pas- 



tur] |1967[ [Sengupta and Mitra[ |1999[ |Tulino and Verdu 
2004). Another way of accounting for these correlations 



is to use the Dyson gas picture in which the possibility for 
two points Yi to be very close to each other can be taken 
into account when obtaining the one-body potential V 9 , 
leading to correlated H ia . 

Another way of extending the analysis of ERMs re- 
viewed here is to use the representation A = HTH^ to 
access quantities that are more advanced than the eigen- 
value density p{A). For Hermitian matrices, the con- 
nected part of the two-point correlation function 

1 N N 

Pc (a, a') = —(Y, S ( A - A «)^ A ' - A »')>- ( 35 °) 



where (xy) c = (xy) — (x)(y), can be obtained from the 



two-point resolvent 
1 

1 



g c (z,z') = 



Tr- 



-A z> - 



A 



(351) 



d z 8 z , (Trlog(z - A)Trlog(*' - A)\ 



(352) 



by invoking the relation (Brezin and Zee 1994): 

p c (A,A') = -^[<7c(+,+) + <?c(- -) 

- <?c(+,-)-Sc(- +)], (353) 

where we introduced a shorthand notation £f c (±j ±) = 
g c (A ± ie, A' ± ie'). Using an elegant diagrammatic ap- 
proach, Brezin and Zee (1995) showed that g c (z,z') can 
be expressed as 

9c(z, z') = -^d z d z , log [1 - U(z, z')g(z)g(z')] , (354) 

where U{z, z') is the irreducible vertex that contains the 
sum of all irreducible diagrams contained in the expan- 
sion of g c (z,z'). For ERMs A = HTH \ the diagram- 
matic method developed in Sec. II. F. 3 can be used to 
obtain 



T 2 



l-g(z)f l-g(z')f 



(355) 



in the limit N — > op, with g{z) the solution of Eq. ( 125 1 . 
Inserting Eq. ( 355 1 into Eq. ( 354 ) yields the two-point 



correlation function ( 353 ) 



Finally, the application of the theory of Euclidean ran- 
dom matrices may be quite fruitful in those branches of 
physics where they naturally appear but were treated 
only numerically until now. In a metal, magnetic mo- 
ments localized at random positions r,j interact indi- 
rectly via polarization of conduction electrons through 
the Ruderman-Kittel-Kasuya-Yosida (RKKY) potential. 
In the three-dimensional space, this interaction can be 
described by an ERM ( |Aristov| |1997[ ) 



J' i 



cos(2fc F |ri - r 3 \) 



(356) 



where kp is the Fermi wave vector. Another example is 
the problem of light scattering in an ensemble of point 
scatterers at positions that one wants to analyze with 
account for the vector character of light. It can be studied 



using an ERM composed of 3 x 3 blocks (Goetschy 2011 
Rusek et al 



1996) 



n 3 exp(ifc ry) 



3i 3 

krij {krij) 2 



1 



1 



(kr r ^- 



(357) 



where ry = r,; 



Tj. We are sure that, besides these two 



examples to which the methods described in the present 
review can be readily applied, other physical problems 
await to be studied in the framework of ERM models. 
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